Virtual target screening to rapidly identify potential protein targets of natural products in drug discovery

1.   Introduction

Mathematical epidemiology constructs mathematical models to understand the dynamics of transmissible diseases, study disease control and elimination, compare strategies of treatments and vaccinations, predict the peak time of an epidemic, and estimate the final size of an outbreak. In the process, mathematical theories on transmission dynamics of infectious diseases are enriched.

It was Kermack and McKendrick who established the celebrated threshold theory to prognosticate the occurrence of an epidemic ^[1,2]. Their threshold was given explicitly in terms of the initial size of the susceptible population (the size of the community). If the initial size of the susceptible population is less than the threshold value, no epidemic can occur; otherwise, it can occur. The use of mass-action law by Kermack and McKendrick was the fundamental for them to address the threshold in terms of the density of the susceptible. Adding vital dynamics from Kermack and McKendrick's epidemic models, Hethcote ^[3,4] showed that the threshold is a global critical value (regardless of the initial conditions) for S-I-S models and S-I-R models (even when vaccination was included). Lajmanovich and York ^[5] and Hethcote and van Ark ^[6] found that the global threshold is valid for multi-group S-I-S models for the spread of gonorrhea, but each group had constant sizes.

Starting in the early 1980's, the threshold for epidemic models has been consistently denoted by $ {{{\cal R}_{0}}} $, which is now formally defined as the expected number of secondary cases produced by a typical infected individual during its entire infectious period in a susceptible population. Though it took many iterations to define $ {{{\cal R}_{0}}} $ ^[7], each step along the way to properly interpret this critical biological number is a prominent landmark in mathematical epidemiology. For instance, basic reproduction rate had been seen in ^[8,9,10]; infectious contact number and replacement number were seen in ^[3,4,11]; the basic reproduction rate later was replaced by basic reproduction ratio such as in ^[12,13]; and finally basic reproduction number has been uniformly accepted since the earlier 1990's ^[14].

As a fundamental threshold, the basic reproduction number plays the foremost role in determining whether a disease can take off or not. If the basic reproduction number for a disease is greater than one, that disease can invade a population. If it's less than one, the disease dies out (for instance see ^[12,13]). The primary purpose of this paper, however, is to demonstrate that the beauty of $ {{{\cal R}_{0}}} $ is not the only story in mathematical epidemiology.

In fact, many epidemiological models have questioned whether a single threshold could completely determine the dynamics of a disease. One has observed that multi-threshold values are necessary to completely determine the dynamics of the disease ^{[15,16,17,18]}. As a result, a collection of epidemic models have shown that bi-stability occurs when $ {{{\cal R}_{0}}} < 1 $, which typically has been driven by the occurrence of backward bifurcations ^{[19,20,21,22]}.

As is shown schematically in Figure 1, when the bifurcation at $ {{{\cal R}_{0}}} = 1 $ is forward, $ {{{\cal R}_{0}}} > 1 $ implies disease can invade population; while $ {{{\cal R}_{0}}} < 1 $, the disease cannot invade the population. If the bifurcation is backward at $ {{{\cal R}_{0}}} = 1 $, there exists a stable equilibrium (even a stable limit cycle exists) when $ {{{\cal R}_{0}}} < 1 $ ^[17], i.e., the disease does not die out even though $ {{{\cal R}_{0}}} < 1 $.

The study of backward bifurcations has a fairly long history. We briefly review the earlier research in mathematical epidemiology on backward bifurcations. Mathematical epidemiology literature first reported the multiple endemic equilibria in mathematical models for the spread of HIV in 1989 by Castillo-Chavez and collaborators ^[23]. It historically challenged that "thinking that S-I-R epidemic models do not have multiple endemic equilibria is not accurate". Their later papers further confirmed that multiple group model of S-I-R type can have multiple endemic equilibria due to the varying sub-population sizes ^[24]. By considering that the birth rate is affected by the disease, Diekmann and Kretzschmar observed multiple stable equilibria in 1991 ^[25]. However, neither ^[23] nor ^[25] explicitly used the words– backward bifurcation. It wasn't until 1995 that Hadeler and Castillo-Chavez first formally used the name: backward bifurcation, and its the occurrence was due to the introduction of core group and non-core group ^[26].

Hadeler and van der Driessche's paper in 1997 entitled "Backward bifurcation in epidemic control" ^[22] studied groups of two that have different susceptibilities can lead to the bi-stability situation. Dushoff et al. explored possible mechanisms behind backward bifurcations in simple disease models in 1998 ^[27]. This research further developed a theoretical local bifurcation indicator that determines whether forward or backward from earlier work in ^[24]. Later, this indicator was well formulated in ^[18] and ^[21]. The new millennium saw an new exploration in the study of backward bifurcations. Variety of mechanisms behind backward bifurcations have been identified in epidemic models. Castillo-Chavez and Feng ^[15] constructed dynamic model for the transmission of tuberculosis, claiming that exogenous reinfections can result in backward bifurcations, which was later found in other mathematical models for the transmission of tuberculosis ^[28,29,30]. It was shown by Kribs-Zaleta and Velasco-Hernández in ^[31] and by Brauer in ^[19,20] that vaccinations can trigger backward bifurcation. Martcheva and Thieme showed that super-infections, a sort of reinfection, are one mechanism in generating backward bifurcations ^[32]. When looking back the mathematical models for tuberculosis ^{[15,28,29,30]} collectively, one can conclude that reinfection would naturally cause backward bifurcations. Wang ^[33] and our previous work ^[17] showed that non-standard treatment strategies can cause the occurrence of backward bifurcation. Here, we very briefly review some well-known studies for backward bifurcations in deterministic differential equation models for disease transmission. It is worth mentioning that backward bifurcations have been found in stochastic and discrete epidemic models (for instance, see ^[34,35]).

The basic reproduction number was found in measuring the primary infection force. When reinfections take place, how do we measure the reinfection forces? Corresponding to $ {{{\cal R}_{0}}} $, this paper finds a basic reinfection number that measures the reinfection forces, which turns out to be a new tipping point for disease dynamics.

The following Theorem 1 is a successful and convenient approach in determining whether or not a system exhibits a backward bifurcation (^[18,21,27]). It will be repeatedly used in this article to guide all our calculations.

Theorem 1. ^[21] Consider a system of ordinary differential equations

with a parameter $ \phi $. Assume that:

1. $ 0 $ is an equilibrium of the system, that is, $ f(0, \phi) \equiv 0 $ for all $ \phi $; and

2. Zero is a simple eigenvalue of $ D_xf(0, 0) = \left(\frac{\partial f_i}{\partial x_j}(0, 0)\right) $ and all other eigenvalues of $ D_xf(0, 0) $ have negative real parts.

Let $ \mathbf{W} = [w_1, w_2, \cdots, w_n]' $ and $ \mathbf{V} = [v_1, v_2, \cdots, v_n] $ be a right and a left eigenvector matrix $ D_xf(0, 0) $, respectively, corresponding to the zero eigenvalue; and let $ f_k(x, \phi) $ be the $ k $th component of $ f(x, \phi) $. Then the local dynamics of system $(1.1)$ around the equilibrium $ 0 $ is totally determined by the signs of $ a $ and $ b $ below:

Particularly, if $ b > 0 $ and $ a > 0 $, then a backward bifurcation occurs for system $(1.1)$ at $ \phi = 0 $; while $ b > 0 $ and $ a < 0 $, a forward bifurcation occurs.

When applying Theorem 1 to an epidemiological model, the equilibrium $ 0 $ corresponds to the disease-free equilibrium and $ \phi = {{{\cal R}_{0}}}-1 $. That is, $ \phi = 0 $ corresponds to $ {{{\cal R}_{0}}} = 1 $. Theorem 1 here is not exactly the same as Thorem 4.1 in ^[21], because we have dropped the assumption that $ \mathbf{W} $ is nonnegative (see Remark 1 in ^[21]).

2.   The Basic Reinfection Number

In this section, we use a simple epidemic model to introduce our essential concept: the basic reinfection number.

We first introduce a mathematical model for a socially transmitted disease such as cigarettes smoking (or drug abuse), where peer pressure or peer influence presumably plays a critical role. The entire population is divided into nonsmokers, smokers, temporary quitters, and permanent quitters. Let $ S $ be the number of nonsmokers; $ I $ the number of smokers; $ R $ the number of the individuals who temporarily quit smoking but have not decided if quit smoking forever; and $ T $ the number of individuals who has commit to quitting cigarettes smoking permanently. A system of ordinary differential equations is used to model the dynamics of cigarettes smoking. The model equations read

The per-capita birth rate and the per-capita death rate are assumed to equal ($ \mu $). $ \beta_1 S \frac{I}{N} $ is the rate of new smokers due to peer pressure. $ \beta_2 R \frac{I}{N} $ is the relapse rate, which is also driven by peer pressure. $ \alpha I $ is the rate of temporary quitting. $ \gamma R $ is the rate of permanent quitting.

In the study of epidemics models, it has become well known that the infection force $ \beta_1 S \frac{I}{N} $ is measured by the basic reproduction number. That is the well-known threshold phenomenon. When $ {{{\cal R}_{0}}} > 1 $, the disease takes off; and $ {{{\cal R}_{0}}} < 1 $ the disease dies out. Now we are interested in how to measure the reinfection force (relapse rate) $ \beta_2 R \frac{I}{N} $ and the role of the reinfection force in determining the dynamics of the disease.

The basic reproduction number for model (2.1) is given by

But the dynamics of the model (2.1) can be described by two quantities: $ {{{\cal R}_{0}}} $ and

Theorem 2. Consider model $(2.1)$. When $ {{{\cal R}_{r}}} < 1 $, a forward bifurcation occurs at $ {{{\cal R}_{0}}} = 1 $; Conversely, when $ {{{\cal R}_{r}}} > 1 $, a backward bifurcation occurs at $ {{{\cal R}_{0}}} = 1 $.

According to Theorem 2, the direction of bifurcation at $ {{{\cal R}_{0}}} = 1 $ is completely determined by $ {{{\cal R}_{r}}} $. We define $ {{{\cal R}_{r}}} = \left(\frac{\beta_2}{\mu+ \alpha}\right)\left(\frac{\alpha}{\mu+ \alpha} \right) \left(\frac{\mu}{\mu+ \gamma} \right) $ to be the basic reinfection number for the SIRT model (2.1). Formally, basic reinfection number, denoted by $ {{{\cal R}_{r}}} $, is defined to be the average number of reinfections of a typical infectious during entire infectious period in a subpopulation of all recovered individuals.

The basic reinfection number $ {{{\cal R}_{r}}} $ in (2.3) has clear epidemiological and biological meaning. Exactly the same as $ {{{\cal R}_{0}}} $, $ {{{\cal R}_{r}}} $ measures both the availability of resources and capability to reinforce. Each term in the expression has clear demographic or epidemiological implication.

Obviously, $ \beta_2 $ is the transmission rate of reinfection; $ \frac{1}{\mu+ \alpha} $ is the average reinfection period which is the same as the average infection period. The first part of the expression in (2.3) is the same as $ {{{\cal R}_{0}}} $, except the for reinfection rate. The second and the third parts imply that the reinfection takes place in the recovered subpopulation ($ R $). To become a member of the $ R $-class, an infected individual must recover from $ I $-class to $ R $-class. This occurs with the probability $ \alpha/(\alpha+\mu) $. Furthermore, an individual in $ R $-class can be reinfected only if the individual does not become a member of $ T $-class. This event happens with the probability $ \mu/(\gamma+\mu) $. Overall $ \left(\alpha/(\alpha+\mu)\right)\left(\mu/(\gamma+\mu)\right) $ gives the probability of contacting the recovered class. Comparing the expressions (2.3) and (2.2), it can be seen that the only difference is that they each apply to distinct sub-populations. The basic reproduction number is applied to the general susceptible population (naive population). The probability of contacting the susceptible class is one because it is assumed that all individuals are susceptible from the definition of $ {{{\cal R}_{0}}} $. However, this is not the case for basic reinfection number, since reinfections take place only if the recovered population is established. The term $ \left(\alpha/(\alpha+\mu)\right)\left(\mu/(\gamma+\mu)\right) $ in (2.3) is exact the probability of the establishment of the recovered population. Therefore, the basic reinfection number is applied to the recovered population. In short, the basic reproduction number takes care of the primary infections and the basic reinfection number accounts for reinfections, which were missing in $ {{{\cal R}_{0}}} $ (It is beyond the reach of $ {{{\cal R}_{0}}} $).

Below, we mathematically prove Theorem 2.

Proof. The Jacobian of system (2.1) around the smoking-free equilibrium $ (N, 0, 0, 0) $ is given by

Choosing $ \beta_1 $ as the bifurcation parameter, then $ {{{\cal R}_{0}}} = 1 $ corresponds to $ \beta_1 = \beta^*_1 = \mu+ \alpha $. It can be seen that $ \lambda = 0 $ is a simple eigenvalue of the matrix $ J_{(\text{DFE,} \beta = \beta^*_1)} $. A right and left eigenvector corresponding to the zero eigenvalue are $ \mathbf{W}' = [w_1, w_2, w_3, w_4] = [- \frac{\beta_1^*}{\mu}, 1, \frac{ \alpha}{\mu+ \gamma}, \frac{ \alpha \gamma}{\mu(\mu+ \gamma)}] $ and $ \mathbf{V} = [v_1, v_2, v_3, v_4] = [0; 1; 0; 0]$, respectively. Using formula of (1.2) and (1.3), we obtain $ b = \frac{\Lambda}{\mu} > 0 $ and $ a = 2\beta_1^*w_1w_2+2 \frac{ \beta_2 \alpha }{\mu + \gamma} $ which can be conveniently expressed as $ \frac{a}{2} = \frac{ \beta_2 \alpha }{\mu + \gamma}- \frac{(\mu+ \alpha)^2}{\mu} $. Therefore, a backward bifurcation occurs if $ \frac{ \beta_2 \alpha }{\mu + \gamma}- \frac{(\mu+ \alpha)^2}{\mu} > 0 $. It is equivalent to $ \left(\frac{\beta_2}{\mu+ \alpha}\right)\left(\frac{ \alpha}{\mu+ \alpha} \right) \left(\frac{\mu}{\mu+ \gamma} \right) > 1 $, i.e. $ {{{\cal R}_{r}}} > 1 $. Applying Theorem 1, Theorem 2 follows.

3.   Examples of Basic Reinfection Number

This section demonstrates the basic reinfection numbers in some epidemic models for tuberculosis and socially transmitted diseases.

3.1.   Exogenous Reinfection of Tuberculosis

For more general tuberculosis (TB) exogenous reinfection models, we refer the reader to ^[15,21]. Here, we construct a simpler model to illustrate the infection number. Basically, we assume that successfully treated individuals are still susceptible individuals. The host population is divided into the following three epidemiological classes: susceptible (S), exposed (L, infected but not infectious), infectious (I), and successfully treated individuals are still susceptible. Let $ N $ denote the total population size. The model takes the following ODE form of three equations:

Here, $ \Lambda $ is the recruitment rate of the population; $ \beta_1 $ is the average number of effective contacts a susceptible has per unit of time; $ \beta_2 $ is the average number of effective contacts an exposed individual has per unit of time. $ \mu $ is the per-capita natural death rate; $ k $ is the primary progression rate at which an exposed individual leaves the latent class by becoming a member of an infectious class; $ d $ is the per-capita disease-induced death rate; $ r_2 $ and $ r_1 $ are the per-capita treatment rates for the exposed and the infected, respectively. The term $ \beta_2L\frac{I}{N} $ models the exogenous reinfection rates. The treated infections and exposed become members of susceptible class. The basic reproduction number for (3.1) is

Because of the exogenous reinfection, model (3.1) undergoes a richer dynamics due to the occurrence of backward bifurcation, which can be characterized by the basic reinfection number.

The basic reinfection number for model (3.1) is

Again, we conclude that the basic reinfection number characterizes the occurrence of backward bifurcation. We formally summarize this result into theorem 3 below. The proof of this theorem can be found in Appendix A.

Theorem 3. Consider model $(3.1)$. The basic reinfection number is given by equation $(3.2)$. When $ {{{\cal R}_{r}}} < 1 $, a forward bifurcation occurs at $ {{{\cal R}_{0}}} = 1 $; Conversely, when $ {{{\cal R}_{r}}} > 1 $, a backward bifurcation occurs at $ {{{\cal R}_{0}}} = 1 $.

Similar to the interpretation to the basic reinfection number for SIRT model (2.1), here is an explanation to the basic infection number for model (3.1). $ \frac{\mu+r_2}{\mu+k+r_2} $ is the fraction of exposed population that does not naturally progress to active TB. This is the target population for exogenous reinfections. Only this portion of exposed population can be reinfected. It was observed that $ \frac{1}{\mu+d+r_1+k} $ tends to be the average reinfection period, which is shorter than the average primary infection period $ \frac{1}{\mu+d+r_1} $.

A perfect match of basic reproduction number and the basic reinfection number can be identified in an extreme case where $ k = 0 $. In this case, $ {{{\cal R}_{0}}} = 0 $ and there is no primary progression and all new cases are the results of the exogenous reinfections. It is amazing that $ {{{\cal R}_{r}}} = \frac{ \beta_2}{\mu+d+r_1} $ has the exactly same formula as the basic reproduction number. If this number is greater than one, then there is possibility that the disease becomes endemic. If the number is less than one, the disease dies out. This study has observed an incredible result: a disease can become endemic even though the basic reproduction number is equal to zero provided that there are reinfections and basic reinfection number is bigger then one.

3.2.   TB Model with Two Strains

Though regular mycobacterium tuberculosis is curable, it takes eight months to clear the bacterium with a complete treatment regime of antibiotics. Incomplete treatment puts the patient at higher risk of developing antibiotic resistance. Now we extend our model in previous subsection to consider multiple strains of tuberculosis infections. An antibiotic sensitive strain of TB and an antibiotic resistant strain are considered. Our goal is to naturally extend our basic reinfection number from one strain to multiple strains.

Notations are the same as model (3.1). Here, we simply add subscripts to distinguish the different strains. The first strain is drug-sensitive and the second is drug-resistant. Let $ L_1 $ and $ L_2 $ denote individuals exposed to regular tuberculosis (antibiotic sensitive) and antibiotic resistant tuberculosis, respectively. Likewise, $ I_1 $ and $ I_2 $ stand for infectious individuals of antibiotic sensitive and antibiotic resistant, respectively. Treatments are applied to both the exposed and the infectious of the sensitive strain. $ p $ is the proportion of exposed individuals of antibiotic sensitive who did not complete their treatment and develop into antibiotic resistance. $ 1-p $ is the proportion of successfully treated exposed individuals. Hence, the treatment rate $ r_2L_1 $ for the exposed class is broken into successful treatment rate $ (1-p)r_2L_1 $ and unsuccessful treatment rate $ pr_2L_1 $. The former becomes member of $ S $ class, the latter goes to $ L_2 $ class. Let $ q $ be the proportion of infectious individuals of antibiotic sensitive who did not complete their treatment and develop into antibiotic resistant; $ 1-q $ be the proportion of successful treatment for infectious individuals. Among total treatment rate of $ r_1I_1 $, $ qr_1 I_1 $ gives the rate at which individuals develop resistant tuberculosis; and $ (1-q)r_1 I_1 $ goes to $ S $ class as the result of successful treatment. The model takes the following form of ordinary differential equations:

Computations in Appendix D show that $ {\cal R}_0^1 = \left(\frac{\beta_1}{\mu+d_1+r_1}\right) \left(\frac{k_1}{\mu+k_1+r_2}\right) $ is the basic reproduction number for drug-sensitive strain and $ {\cal R}_0^2 = \left(\frac{\beta_1}{\mu+d_2}\right)\left(\frac{k_2}{\mu+k_2}\right) $ is the basic reproduction number for drug-resistant strain. Then the basic reproduction number for whole system is

In Appendix D, a detailed computation for formula (3.4) is presented.

The situation is similar to that of the basic reproduction number (3.4). For each strain of TB, there is a corresponding basic reinfection number. The basic reinfection number for sensitive strain is

which is derived from Theorem 1 by assuming $ {\cal R}_0^1 > {\cal R}_0^2 $. On the other hand, if $ {\cal R}_0^1 < {\cal R}_0^2 $, the basic reinfection number for resistant strain of TB is

Hence, the basic reinfection number for the entire system is

Theorem 4. Consider model $(3.3)$. The basic reinfection number is given by $(3.5)$. When $ {{{\cal R}_{r}}} < 1 $, a forward bifurcation occurs at $ {{{\cal R}_{0}}} = 1 $; Conversely, when $ {{{\cal R}_{r}}} > 1 $, a backward bifurcation occurs at $ {{{\cal R}_{0}}} = 1 $.

We conclude again that the basic reinfection number controls the occurrence of backward bifurcation. Biological interpretation for the basic infection number is exactly the same as (3.2) of TB model (3.1) for a single strain. Both resistant strain and sensitive strain include the exogenous reinfection which causes the occurrence of backward bifurcations. What if both strains generate bi-stability? Numerically, tri-stability dynamics were found in ^[28]. That is, a collision of two backward bifurcations is able to generate tri-stability dynamics. We refer the interested readers to the examples in ^[28].

3.3.   A Mathematical Model for Sleeper Effects

Children who smoked once before age 11 are twice as likely to become a regular smoker by age 14 ^[37]. This is so called the sleeper effects. In this subsection, we present a mathematical model to study the dynamics of smoking in the United States among children ages 11 to 18 ^[38].

All children enter the system as susceptible non-smokers with either little or no smoking history. Consequently, a two-string SIRT model is constructed, where string $ 1 $ is composed of individuals who smoked once before age 11, and string $ 2 $ is composed of those who did not. Using this model, the "sleeper effect" can be investigated.

Each string is composed of four classes; $ N_i $, $ S_i $, $ Q_i $ (where $ i = 1, 2 $ denotes the respective string), and $ P $. The susceptible non-smoker classes ($ N_i $) are composed of individuals who have never been regular smokers. The regular smoker classes ($ S_i $) are composed of smokers. The recovered classes ($ Q_i $) are composed of individuals who have quit smoking but are susceptible to relapse. The class of permanent non-smokers ($ P $) is composed of individuals who have decided to permanently quit (from $ Q_i $).

Individuals who turn 11 years old enter either $ N_1 $ or $ N_2 $ with rates $ \mu rT $ and $ \mu \left(1-r \right)T $, respectively, where $ T $ is the total population of adolescents. Non-smokers ($ N_i $) are susceptible to becoming regular smokers as a result of interactions with regular smokers ($ S_i $). Smokers ($ S_i $) can quit smoking due to personal choice or education. Individuals who have quit ($ Q_i $) can relapse due to social interactions with smokers. All individuals who are not currently regular smokers ($ N_i $, $ Q_i $) can move to class $ P $ by choosing to permanently abstain from cigarettes. We assume the total population, $ T $, is constant. Definitions of all model parameters are shown in Table 1. The following mathematical model was proposed in ^[38].

The equations of $ N_1 $, $ S_1 $ and $ Q_1 $ in model (3.6) describe the dynamics of the population that smoked at least once before age 11. Likewise, the equations of $ N_2 $, $ S_2 $ and $ Q_2 $ describe the dynamics of the population that did not smoke before age 11.

The basic reproduction number for model (3.6) is

$ \frac{\beta_1}{\phi_1+\mu} $ is the basic reproduction number of subpopulation with a smoking history and $ \frac{\beta_2}{\phi_2+\mu} $ is the basic reproduction numbers of subpopulation without a smoking history. The overall basic reproduction number is the weighted average of the basic reproduction numbers for these two populations using $ r $ and $ (1-r) $ as the weights.

The general structure of model (3.6) is almost the same as model (2.1). Since each string of (3.6) has the same structure as (2.1), we then expect to observe the same structure of basic reinfection number. To highlight the role of the basic reinfection number in this research, we here consider a simpler case where that $ \beta_1 = \beta_2 $.

The basic reinfection number for model (3.6) is

where $ \left(\frac{\psi_1}{\phi_1+\mu}\right)\left(\frac{\phi_1}{\phi_1+\mu}\right) \left(\frac{\mu}{ \alpha_1+\mu}\right) $ is the basic reinfection number of subpopulation with a smoking history before age 11. Likewise, $ \left(\frac{\psi_2}{\phi_2+\mu}\right) \left(\frac{\phi_2}{ \alpha_2+\mu}\right)\left(\frac{\mu}{ \alpha_2+\mu}\right) $ is the basic reinfection number of subpopulation without a smoking history before age 11.

To interpret the meaning of $ \frac{\psi_1}{\phi_1+\mu} \frac{\phi_1}{\phi_1+\mu} \frac{\mu}{ \alpha_1+\mu} $, we break it into two components. $ \frac{\psi_1}{\phi_1+\mu} $ is the average of relapse due to a typical smoker with smoking history before age 11 since $ \psi_1 $ is the relapse rate and $ \frac{1}{\phi_1+\mu} $ is the average smoking period in $ S_1 $. The second part $ \frac{\phi_1}{\phi_1+\mu} \frac{\mu}{ \alpha_1+\mu} $ finds the proportion of potential relapse population, because $ \frac{\phi}{\phi_1+\mu} $ is the proportion of quitting in $ S_1 $, and $ \frac{\mu}{\alpha_1+\mu} $ is the proportion of not permanent abstaining from cigarettes smoking of $ Q_1 $ class.

The overall basic reinfection number (3.8) is the weighted average of the basic reinfection numbers of two subpopulations using $ r $ and $ (1-r) $ as the weights. It has the same format as the basic reproduction number (3.7).

We have already achieved the same conclusion that if $ {{{\cal R}_{r}}} > 1 $, the model (3.6) undergoes a backward bifurcation at $ {{{\cal R}_{0}}} = 1 $. Such a statement is justified rigorously in Appendix C.

If relapse does not come into play, i.e., $ \psi_i = 0 $, then $ {{{\cal R}_{r}}} = 0 < 1 $ and consequently backward bifurcation cannot happen. Namely, $ {{{\cal R}_{r}}} > 1 $ happens uniquely because of the relapse. To handle with cigarette smoking in adolescence, relapse must be avoided by all means.

4.   Conclusions

The reinfection forces (or capability of relapse) is successfully measured by the basic reinfection numbers, the identical role that basic reproduction number plays for the primary infection forces. A primary infection invades a population if the basic reproduction number is greater than one. If reinfection or relapse exist, the basic reinfection indicates a successful invasion of a disease into a population even when the basic reproduction number is less than one. Not only does the basic reinfection number $ {{{\cal R}_{r}}} $ quantify the reinfection forces, but also specifies the type of bifurcation at the basic reproduction number $ {{{\cal R}_{0}}} = 1 $.

In reconsidering the remarkable observation from model (3.1) for tuberculosis, where the disease is still able to spread when the basic reproduction number is equal to zero, we may want to combine $ {{{\cal R}_{0}}} $ and $ {{{\cal R}_{r}}} $ together for a complete elimination of the disease. For instance, max$ \{{{{\cal R}_{0}}}, {{{\cal R}_{r}}}\} $ would be a more proper control number, since the disease dies out when max$ \{{{{\cal R}_{0}}}, {{{\cal R}_{r}}}\} < 1 $ regardless the initial epidemiological status. On one hand, max$ \{{{{\cal R}_{0}}}, {{{\cal R}_{r}}}\} $ can be applied to epidemic models without reinfections because, when there is no reinfection, the basic reinfection number $ {{{\cal R}_{r}}} = 0 $, subsequently max$ \{{{{\cal R}_{0}}}, {{{\cal R}_{r}}}\} = {{{\cal R}_{0}}} $. On the other hand, more important observation is that max$ \{{{{\cal R}_{0}}}, {{{\cal R}_{r}}}\} $ is able to explain why an infectious disease can spread when $ {{{\cal R}_{0}}} = 0 $. Therefore, as an indicator whether the disease can take off or not, max$ \{{{{\cal R}_{0}}}, {{{\cal R}_{r}}}\} $ is better than $ {{{\cal R}_{0}}} $ alone.

Using a few of deterministic epidemic models, we intend to introduce a new concept in theoretic epidemiology. The basic reinfection number $ {{{\cal R}_{r}}} $ has successfully characterized the occurrence of backward bifurcations for a large number of epidemic models. But we found that it is unsuccessful for the appearance of backward bifurcations caused by vaccination factor or nonlinear treatment rate^[17,31,33]. Further research would find a corresponding counterpart for the epidemiological model where vaccination matters or nonlinear treatments are applied.

Unlike $ {{{\cal R}_{0}}} $, which has been well formulated by using the next generation operator ^[13,18,36], in this paper we do not have a universal formula for $ {{{\cal R}_{r}}} $. Our approach to formulate it is to compute $ a $ and $ b $ from Theorem 1. Then manage the inequality $ a > 0 $ in order to find the basic reinfection number. Future work may be devoted to formulate $ {{{\cal R}_{r}}} $ in general scenarios.

In the range of $ {{{\cal R}_{0}}} < 1 $, if reinfection force is strong enough to make $ {{{\cal R}_{r}}} > 1 $, the disease may be persistent even when the primary infection cannot support an endemic. This does not work for the full range of $ {{{\cal R}_{0}}} < 1 $, because when the basic reproduction number is too small, there are not enough recovered individuals to be reinfected. That is why any backward bifurcation curve eventually has to turn around, as can be seen in Figure 4.

Conflict of interest

The author declares that there is no conflict of interest regarding the publication of this paper.

Appendix

A.   Proof of Theorem 3

Here we show how to get the reinfection number for model (3.1). The Jacobian of the system (3.1) around the TB-free equilibrium is

When $ {{{\cal R}_{0}}} = \left(\beta_1/(\mu+d+r_1)\right)\left(k/(\mu+k+r_2)\right) = 1 $, a left eigenvector of $ J $ is $ \mathbf{V} = \left[

$$\begin{matrix} v_1 & v_2 & v_3\end{matrix}$$ \right] = \left[ $$\begin{matrix} 0 & v_2 & v_2 (\mu+k+r_2)/k\end{matrix}$$ \right] $ and a right eigenvector is $ \mathbf{W}' = \left[w_1, w_2, w_3\right] $ $ = \left[w_2(r_2+(\mu+k+r_2)(r_1-\beta_1)/\beta_1)/\mu, w_2, w_2(\mu+k+r_2)/\beta_1 \right] $. Choosing $ v_2 = w_2+1 $ with any $ w_2 > 0 $, we obtain that

To single out the sign of $ a $, we rewrite the last equation as

Obviously, $ a > 0 $ if and only if $ \left(\frac{ \beta_2}{\mu+d+r_1+k}\right)\left(\frac{\mu+r_2}{\mu+k+r_2}\right) > 1 $. Then

is the basic reinfection number. Applying Theorem 1, we conclude our proof.

B.   Computation of $ {{{\cal R}_{r}}} $ for Model (3.3)

We choose $ \beta_1 $ as the bifurcation parameter. Denote the right-hand functions of (3.3) by $ f_1 $, $ f_2 $, $ f_3 $, $ f_4 $, $ f_5 $; and the state variables $ (S, L_1, L_2, I_1, I_2) $ are denoted by $ (x_1, x_2, x_3, x_4, x_5) $. Set $ K = \frac{\Lambda}{\mu} $ for a shortcut.

Applying theorem 1, we directly compute $ a $ and $ b $. The Jacobian of the system evaluated at the TB-free equilibrium $ (\Lambda/\mu, 0, 0, 0, 0) $ is

Comparing with one-strain TB model (3.1), we here have to additionally deal with whether $ {{{\cal R}_{0}}} = {\cal R}_0^1 = \left(\frac{\beta_1}{\mu+d_1+r_1}\right) \left(\frac{k_1}{\mu+k_1+r_2}\right) $ or $ {{{\cal R}_{0}}} = {\cal R}_0^2 = \left(\frac{\beta_1}{\mu+d_2}\right)\left(\frac{k_2}{\mu+k_2}\right) $ because the overall basic reproduction number is $ {{{\cal R}_{0}}} = \text{max} \left\{{\cal R}_0^1, {\cal R}_0^2 \right\} $. For convenience, we tabulate the non-zero second derivatives ($ \partial^2f_k/\partial x_i \partial x_j $, $ i, j, k = 1\cdots 5 $) evaluated at the TB-free equilibrium $ (\Lambda/\mu, 0, 0, 0, 0) $ and $ {{{\cal R}_{0}}} = 1 $ in Table 2.

If $ {\cal R}_0^1 > {\cal R}_0^2 $, then $ {{{\cal R}_{0}}} = 1 $ means $ {\cal R}_0^1 = 1 $. A left eigenvector of (B.1) at $ {\cal R}_0^1 = 1 $ is $ \mathbf{V} = [v_1, v_2, v_3, v_4, v_5] = [0, \frac{k_1}{\mu+r_2+k_1}v_4, 0, v_4, 0] $, where $ v_4 $ is free. A right eigenvector is $ \mathbf{W} = \left[w_1, w_2, w_3, w_4, w_5\right]' $, where

It can be checked $ b = \frac{k_1}{\mu+r_2+k_1}v_4w_4 > 0 $. To compute $ a $, we need to pay attention to $ f_2 $ and $ f_4 $ only because $ v_1 = v_3 = v_5 = 0 $.

Looking up Table 2, we obtain that

Therefore, $ a > 0 $ iff $ \beta_2(\mu+r_2) > (\mu+d_1+r_1+k_1)(\mu+k_1+r_2) $, i.e. $ \frac{\beta_2}{\mu+r_1+d_1+k_1} \frac{\mu+r_2}{\mu+k_1+r_2} > 1 $. Then the basic reinfection number for sensitive strain is

Now assume $ {\cal R}_0^1 < {\cal R}_0^2 $. Then $ {{{\cal R}_{0}}} = 1 $ means $ {\cal R}_0^2 = 1 $. A left eigenvector of (12.1) at $ {\cal R}_0^2 = 1 $ is $ \mathbf{V} = [v_1, v_2, v_3, v_4, v_5] $, where $ v_1 = 0 $, $ v_2 = \frac{pr_2(\mu+d_1+r_1)+k_1qr_1}{k_2(\mu+d_1+r_1)(\mu+k_1+r_2)-k_1(\mu+d_2)(\mu+k_2)} \frac{k^2_2}{\mu+d_2}v_5 $, $ v_3 = \frac{k_2}{\mu+d_2}v_5 $, $ v_4 = \frac{k_2qr_1(\mu+k_1+r_2)+pr_2(\mu+d_2)(\mu+k_2)}{k_2(\mu+d_1+r_1)(\mu+k_1+r_2)-k_1(\mu+d_2)(\mu+k_2)} \frac{k_2}{\mu+d_2}v_5 $, and $ v_5 $ is free. A right eigenvector is $ \mathbf{W}' = \left[w_1, w_2, w_3, w_4, w_5\right] = \left[

$$\begin{array}{ccccc} \frac{(\mu+k_2)(\mu+d_2)}{\mu k_2}w_5, & 0, & \frac{\mu+d_2}{k_2}w_5, & 0, & w_5 \end{array}$$ \right] $. $ w_5 $ is free. $ w_2 = w_4 = 0 $ makes the computation of $ a $ a lot easier.

Therefore, $ a > 0 $ iff $ \mu\beta_2 > (\mu+k_2)(\mu+d_2+k_2) $, i.e. $ \frac{\beta_2}{\mu+d_2+k_2} \frac{\mu}{\mu+k_2} > 1 $. Then the basic reinfection number for resistant strain is

The basic reinfection number for entire system is

C.   Computation of $ {{{\cal R}_{r}}} $ for Model (3.6)

The basic reinfection number for model (3.6) is derived in a simple case $ \beta_1 = \beta_2 = \beta $. We look at the Jacobian of the system evaluated at the smoker-free equilibrium $ (rT, (1-r)T, 0, 0, 0, 0, 0) $

Choose $ \beta $ as the bifurcation parameter. Apparently $ {{{\cal R}_{0}}} = 1 $ corresponds to $ \beta = \beta^* = \frac{(\phi_1 + \mu)(\phi_2 + \mu)}{r(\phi_2 + \mu)+(1-r)(\phi_1 + \mu)} $. A right eigenvector associated with the eigenvalue $ \lambda = 0 $ of $ J $ is $ \mathbf{W}' = \left[w_1, w_2, w_3, w_4, w_5, w_6, w_7 \right] $, where

A left eigenvector associated with the eigenvalue $ \lambda = 0 $ of $ J $ is $ \mathbf{V} = [v_1, v_2, v_3, v_4, v_5, v_6, v_7] = \left[0, 0, \frac{\phi_2+\mu}{\phi_1+\mu}, 1, 0, 0, 0\right] $. Without confusion, we set $ f_i $ to be the right-hand side functions of system (3.6) and denote $ [N_1, N_2, S_1, S_2, Q_1, Q_2, P] $ by $ [x_1, x_2, x_3, x_4, x_5, x_6, x_7] $. Because of $ v_1 = v_2 = v_5 = v_6 = v_7 = 0 $, only $ f_3 $ and $ f_4 $ matter when calculating $ a $ and $ b $.

Applying Theorem 1 to model (3.6), we obtain that

Organize the above expression to check the sign of $ a $. Noticing $ w_3 > 0 $ and $ w_4 > 0 $.

It can be seen that $ a > 0 $ if and only if

which is equivalent to

Consequently, the basic reinfection number for model (3.6) is

D.   Computation of $ {{{\cal R}_{0}}} $ for Model (3.3)

The approach of next-generation operator ^[18] is used to calculate the basic reproduction number $ {{{\cal R}_{0}}} $ for model (3.3).

$ X = [L_1, L_2, I_1, I_2] $, $ \frac{dX}{dt} = F(X)-V(X) $, where

$ F(X) = \left[

$$\begin{matrix} \frac{\beta_1 I_1K}{N}\\ \frac{\beta_2 I_2K}{N}\\ 0\\ 0\\\end{matrix}$$ \right] $ and $ V(X) = \left[ $$\begin{matrix} \frac{\beta_2I_1}{N}L_1+(\mu+k_1+r_2)L_1\\ -pr_2L_1-qr_1I_1+\frac{\beta_2I_2}{N}L_2+(\mu+k_2)L_2\\ -\frac{\beta_2I_1}{N}L_1-k_1L_1+(\mu +d_1+r_1)I_1\\ -\frac{\beta_2I_2}{N}L_2-k_2L_2 +(\mu+d_2)I_2\\ \end{matrix}$$ \right]. $

$ {\cal F} = \left. \frac{\partial F}{\partial X}\right|_{X = 0} = \left[

$$\begin{matrix} 0 & \beta\\0 & 0\end{matrix}$$ \right] $ with $ \beta = \left[ $$\begin{matrix} \beta_1 & 0\\0 & \beta_1\end{matrix}$$ \right] $, and $ {\cal V} = \left.\frac{\partial V}{\partial X}\right|_{X = 0} = \left[ $$\begin{matrix} v_{11} &v_{12}\\ v_{21}&v_{22}\end{matrix}$$ \right] $, where $ v_{11} = \left[ $$\begin{matrix} \mu+k_1+r_2 & 0\\ -pr_2 & \mu+k_2\end{matrix}$$ \right] $, $ v_{12} = \left[ $$\begin{matrix} 0 & 0\\-qr_1 & 0\end{matrix}$$ \right] $, $ v_{21} = \left[ $$\begin{matrix} -k_1 & 0\\ 0&k_2\end{matrix}$$ \right] $, and $ v_{22} = \left[ $$\begin{matrix} \mu+d_1+r_1 & 0\\ 0 & \mu+d_2\end{matrix}$$ \right] $. A block format of $ {\cal V}^{-1} = \left[ $$\begin{matrix} \gamma_{11} & \gamma_{12}\\ \gamma_{21} & \gamma_{22}\end{matrix}$$ \right] $. Since $ {\cal F}{\cal V}^{-1} = \left[ $$\begin{matrix} \beta\gamma_{21} & \beta\gamma_{22}\\ 0 & 0\end{matrix}$$ \right] $, the spectrum radius of $ {\cal F}{\cal V}^{-1} $ is determined by $ \beta\gamma_{21} $. Using block operations of matrices, one can find

$ \gamma_{21} = \left[

$$\begin{matrix} \frac{k_1}{(\mu+d_1+r_1)(\mu+k_1+r_2)} & 0\\ \frac{k_1\left(pr_2+ \frac{k_1qr_1}{\mu+d_1+r_1}\right)}{(\mu+d_1+r_1)(\mu+k_1+r_2)(\mu+k_2)} & \frac{k_2}{(\mu+d_2)(\mu+k_2)}\end{matrix}$$ \right] $. $ \beta\gamma_{21} = \left[ $$\begin{matrix} \frac{\beta_1k_1}{(\mu+d_1+r_1)(\mu+k_1+r_2)} & 0\\ \frac{\beta_1k_1\left(pr_2+ \frac{k_1qr_1}{\mu+d_1+r_1}\right)}{(\mu+d_1+r_1)(\mu+k_1+r_2)(\mu+k_2)} & \frac{\beta_1k_2}{(\mu+d_2)(\mu+k_2)}\end{matrix}$$ \right] $. Therefore

and formula (3.4) is proved.