A new perspective on infection forces with demonstration by a DDE infectious disease model

1.   Introduction

Incidence rate in an infectious disease model plays a crucial role in the disease dynamics. It is the rate (i.e., number of cases per unit time) of incidence of new infections coming from the susceptible population. Expressing the incidence rate in the form $f S$ with $S$ being the susceptible population, then the role of the incidence rate is represented by $f$, which is often referred to as the infection force. An infection force can be dependent on $I(t)$ or on both $I(t)$ and $S(t)$ where $I(t)$ and $S(t)$ are the populations of the infectious and susceptible classes of the host. Actually, there have been many works in literature that investigate the disease dynamics described by various types of ODE disease models (SI, SIS, SIR, SIRS, etc.) adopting various infection force functions (or incidence rates), see, e.g., ^{[1,2,3,4,5,6,7,8,9,10,11,12,13]} and the references therein. For example, in Liu et al. ^[6,7] the incidence rate $\beta I^pS^q$ were used for some in SIRS and SEIRS models respectively, and it was later also used by Korobeinikov et al. ^[4]. This incidence rate corresponds to an infection force function of the form $f = f(I, S) = \beta I^p S^{q-1}$. In Ruan et al. ^[10], a saturated infection force $f = f(I) = aI^2/(b+I^2)$ was adopted in an SIRS model. In Korobeinikov et al. ^[5], general incidence rate of the form $g(S, I, N)$, where $N$ is the total population, was discussed and the global properties of some infectious disease models with various nonlinear $g$ were explored. In Wang ^[11], some general nonlinear infection forces were discussed for an SIRS ODE model with balanced demography, and these can include the form $f = \lambda I/(1+pI+qI^2)$ and $f = \lambda I^2/(1+pI+qI^2)$ which are infection forces discussed in, e.g., ^[8,9,12]. In Liu et al. ^[14] for an SEIHR ODE model, an infection force of the form $\beta I e^{-a_1E-a_2I-a_3H}$ was used where $E(t)$ and $H(t)$ are the numbers of exposed and hospitalized hosts; and along the same line but for an SEI ODE model with logistic demography for the host, an infection force of the form $\beta I e^{-m I}$ was used in reference [15]. Particularly in McCallum et al. ^[13], some valuable discussions on transmission functions are presented from biological perspectives, and seven particular transmission functions were collected and compared, including some of the above mentioned ones.

It has been found that some nonlinear infection force functions, particularly those non-monotone ones, can lead to very complicated disease dynamics, as opposed to the linear infection force arising from mass action incidence rate. For details in this regard, a reader is referred to, e.g., ^{[4,8,9,10,12,14,15]} and the references therein.

We note that aforementioned literatures on nonlinear infection force mainly focus on pure theoretical research in dynamical systems, and hence, have brief justifications from the viewpoint of infectivity, and some are even ambiguous and out of reasonable biological explanations. Here we firstly present a novel alternative perspective for considering infection force functions which can not only reconstruct the infection force function mentioned above, but also help motivate our new model. The main idea is that during an epidemic or pandemic, particularly for a newly emerged disease such as the SARS in 2003 and Covid-19 broken out in $2020$, due to various control measures including lockdowns, quarantines and isolations, ordered stay-home and even the effective use of personal protection equipments (PPEs), only a proportion of the susceptible individuals will actually be possibly exposed to the infectious hosts and hence be possibly infected. For convenience of later references, let us call this proportion of susceptible individuals practically susceptible individuals. A susceptible host who takes extreme precaution by either staying home or using full PPEs has zero or little chance to be infected. The practically susceptible population can be expressed as $S_{p} = P S(t)$, where $P \in [0, 1]$ is a measurement of the proportion; alternatively, this $P$ can also be explained as the probability that a susceptible individual is a practically susceptible host. In general, such a proportion is non-pharmaceutical and non-biological, and it depends on the level of precaution that is impacted by the level of severity of the epidemics, mediated by, e.g., the control measures implemented by the government, regulations and guidelines imposed by the public health agency, as well as media coverage. Let $L(t)$ denote the measurement of the level of severity, and accordingly $P = P(L)$. As the severity level increases, people will generally become more precious and more protective by either restricting their activities or using PPEs; moreover, the government may implement more and more restrictive control measures, which will make fewer and fewer people available for infection. These facts justify the following basic assumption for the proportion function $P(L)$:

(H1) $P(L)$ is a monotonic non-increasing function, satisfying $P(0) \leq 1$ and $\lim_{L \rightarrow \infty} P(L) \geq 0$.

Below are just three prototypes of such a function satisfying (H1):

(A) $P_{1}(L) = \frac{m_{1} L+1}{m_{2} L+1}$ where $m_{1}, m_{2}$ are all positive constants satisfying $\frac{m_{1}}{m_{2}} < 1$;

(B) $P_{2}(L) = \frac{m_{1} L+b}{c L^{2}+m_{2} L+b}$, where all parameters are positive constants and satisfy $\frac{m_{1}}{m_{2}} < 1$;

(C) $P_{3}(L) = e^{-h L}$, where $h > 0$.

Note that although $P_j(L), \, \, j = 1, 2, 3$ all satisfy (H1), their decay rates are different, with $P_1(L)$ decaying slowest and $P_3(L)$ decaying the fastest. Also, $P_1(L)$ has a positive lower bound $m_1/m_2 < 1$, and this may accommodate the scenario of maintaining essential services open during a severe epidemic, such as Covid-19.

Now let us consider the mass action infection mechanism which has the incidence rate at time $t$ as $\beta I(t) S(t)$. This corresponds to a linear infection force $f(t) = \beta I(t)$, where $\beta$ is the transmission rate. With the above observation of practically susceptible population, the incidence rate under the mass action infection mechanism should be revised to $\beta I(t) S_p(t) = \beta I(t) [P(L(t))S(t)] = [\beta I(t) P(L(t)) ]S(t)$, resulting in the infection force $f(t) = \beta I(t)P(L(t))$. This way, the infection force at time $t$ is explicitly related to the proportion function $P(L(t))$ which depends on the severity level $L(t)$ at time $t$. When $L(t)$ is measured just by the number of infected individuals at time $t$, i.e., $L(t) = I(t)$, those infection force functions used in the aforementioned literatures can be easily obtained by properly chosen proportion function $P(I)$ for the practically susceptible population. In this case, the infection force $f = f(I) = I P(I)$ can be monotone, either increasing and unbounded, or increasing and saturated if $P(I)$ decrease slowly; it can also be non-monotone (one hump shape) if $P(I)$ decreases fast enough, including those used in references [8,9,10,12,14,15] where complicated dynamics have been observed. We point out that Liu et al. ^[14] proposed and analyzed an SEIR-H model with $E(t)$ and $H(t)$ representing the populations of exposed and hospitalized groups, and the infection force there is a result of adopting $L(t) = a_1 E(t)+a_2 I(t)+ a_3 H(t)$ and $P(L) = e^{-L(t)}$.

Note that in practice, a more reasonable and meaningful measurement for the severity of an epidemic should contain information not only at the present time but also certain information over a period of the past time. When determining their precaution and protection levels, people do look at the numbers of infected cases both at the present time and in recent times. This suggests the following form for $L(t)$:

where the constant $\tau > 0$ specifies a length of time interval that one wants to look at and the weight function $w(\xi)$ reflects the variation of the impact of disease surveillance in the past interval $[t-\tau, t]$ on the severity at the present time. Considering the fact that in reality, the data are typically collected at discrete times (e.g., hourly, daily and weekly), $w(\xi)$ is accordingly taken as $w (\xi) = \sum_{j = 0}^{n}k_j\delta(\tau_j)$, leading to the following discrete form for $L(t)$:

where $k_j \ge 0, \, \, j = 0, 1, \cdots, n$ are the discrete weights, satisfying the normal condition $\sum_{j = 0}^n k_j = 1$; $\tau_0 = 0$ and $0 < \tau_1 < \tau_2 < \cdot < \tau_n$ are the positive numbers that may account for the past $n$ time points at which the infected cases are reported/released. When the data in these past $n$ time points and the present time are all given the same weight, there holds $k_j = 1/(n+1), \, \, j = 0, 1, \cdots, n$.

To demonstrate the above new perspective for infection force, we incorporate the above new framework for infection force into the classic SIRS model to obtain the following model system for the transmission dynamics of an infectious disease of SIRS type:

where, as discussed above, the infection force at time $t$ is given by

with the proportion function satisfying (H1). Although $P(L)$ does not have to be smooth and even continuous to reflect the situation when the control measures are often implemented in terms of time intervals, for convenience of theoretical analysis, we assume $P(L)$ is differentiable, in addition to (H1). Notice that all parameters are positive in model (1.3); in addition to the transmission rate $\beta$, the others respectively stand for

$-$    $\Lambda$: the recruitment rate of the population;

$-$    $d$: the natural death rate of the population;

$-$   $r$: the recovery rate of infective individuals;

$-$    $\epsilon$: the disease-induced death rate;

$-$    $\alpha$: the rate of removed individuals who lose immunity and return to susceptible class.

We point out that in a very special case: $\epsilon = 0, \, \, \alpha = 0$ and $n = 1, \, \, k_0 = 0$ and $k_1 = 1$ (i.e., severity at time $t$ is measured by $L(t) = I(t-\tau)$) and $P(L) = e^{-h L}$, system (1.3) reduces to the delayed SIR model

which was thoroughly discussed recently by Song et al. ^[16], including a local and global bifurcation analysis by using the delay $\tau$ as the bifurcation parameter, which reveals onset and termination of Hopf bifurcation from a positive equilibrium. The authors attribute the term $e^{-h I(t-\tau)}$ as the impact of media.

In more recent work, in further studying the impact of media on the disease dynamics of an SEIS model, Song et al. ^[17] introduced a variable $M(t)$ to denote the average number of news items related to the disease outbreak, and proposed the following SEIS-M model

Threshold dynamics and bifurcation from an endemic equilibrium were explored in reference [17]. Here the media coverage variable $M(t)$ plays a similar role to that of the severity variable $L(t)$ in our model (1.3).

The models (1.3)-(1.4) and (1.5) have something in common: the infection forces contain time delay. We would like to draw readers' attention to another scenario for delay$-$dependent infection forces but for vector-borne diseases. The key idea traces back to Cooke ^[18] where it was assumed that the population of the infectious vectors at time $t$ is proportional to the population of the infectious host in the previous time $t-\tau$ with $\tau$ being the latency of disease in the vector. This suggests the incidence rate for the host $F(I_v(t)) S_h(t) = F(p I_h(t-\tau)) S_h(t)$, and accordingly, reduces a model system for vector-borne disease to a system only containing the respective variables (susceptible and infectious, etc.) for the host populations, but with delayed infection force function. Along this time, there have been many works with the corresponding models typically demonstrating the global threshold dynamics between the disease free equilibrium and endemic equilibrium, thanks to the powerful Lyapunov method. The success of this method is mainly attributed to assuming that $F$ only depends on the delayed term $I_h(t-\tau)$ and $F(\cdot)$ is nondecreasing. For some details, see e.g., ^{[19,20,21,22,23,24,25]} and the references therein.

In the rest of this paper, we explore the dynamics of the SIRS model (1.3) with the infection force given by (1.4). In Section 2, we address the well-posedness of model (1.3)-(1.4) with general $P(L)$. In Section 3, we identify the basic reproduction number $R_0$ and discuss the equilibria of the model system in the desired region. In Section 4, we consider a particular $P(L(t)) = e^{-hL}$ which leads to the special infection force $f(t) = \beta I(t)e^{-h(k_0I(t)+k_1I(t-\tau))}$, resulting to the model system (4.1). For this infection force, and in the case of $R_0 > 1$ so that a unique endemic equilibrium $E^*$ exists, we analyze the stability of $E^*$ and obtain conditions for Hopf bifurcation to occur. In Section 5, we compare our results with those in reference [16] which is a special case of (4.1) ($\alpha = 0, \epsilon = 0, k_0 = 0, k_1 = 1$). Such a comparison particularly reveals the difference of the disease dynamics between using multiple time disease surveillance (i.e., $I(t)$ and $I(t-\tau)$) and only using single time disease surveillance (i.e., $I(t-\tau)$ only). In Section 6, we present some numerical simulations to illustrate the theoretical results. Finally, we conclude the paper by a brief summary together with some discussion, in Section 7.

2.   Well-posedness of the model

As is customary, we denote by $C = C\left([-\tau, 0], \mathcal{R}^{3}\right)$ the set of all continuous functions defined on $[-\tau, 0]$, and by $C_+ = C\left([-\tau, 0], \mathcal{R}_+^{3}\right)$ the subset of all of non-negative functions in $C$. Then $C$ is a Banach space with the maximum norm and $C_+$ is the positive cone of $C$ which induces the natural order in $C$.

By the fundamental theory of functional differential equations (see, e.g., Hale et al.^[26]), for $\phi = \left(\phi_{1}, \phi_{2}, \phi_{3}\right)^{T} \in C$, the model system (1.3) has a unique solution $(S(t), I(t), R(t))$ satisfying

which exists in a maximal interval of existence $(0, T_m)$.

By the nature of the variables in the model (1.3)-(1.4), we are only interested in non-negative solutions. Accordingly, we further require the initial functions to be in $C_+$. For such non-negative initial functions, we show below that the corresponding solution is also non-negative. Firstly, from the second equation in model (1.3)-(1.4), we have

Similarly, the third equation in (1.3) lead to

and thus, by $R(0) \ge 0$ and the non-negativity of $I(t)$ confirmed above, we conclude $R(t) \ge 0$ for $t \in (0, T_m)$. To verify that $S(t) > 0$ for $t > 0$, we assume the contrary. Let $t_{1} > 0$ be the first time such that $S\left(t_{1}\right) = 0,$ then by the first equation of (1.3) we have $S^{\prime}\left(t_{1}\right) = \Lambda+\alpha R(t_1) > 0$ and therefore, $S(t) < 0$ for $t \in\left(t_{1}-\epsilon_1, t_1\right)$ where $\epsilon_1 > 0$ is sufficiently small. This contradicts $S(t) > 0$ for $t \in\left[0, t_{1}\right)$, and the contradiction proves $S(t) > 0$ for all $t \in (0, T_m)$.

Next we prove the boundedness of the solution. To this end, we consider $N(t) = S(t)+I(t)+R(t)$. Simple calculations lead to

This implies that $\lim \sup _{t \rightarrow \infty} N(t) \leq \Lambda/d$, concluding the boundedness of $N(t)$. This together with the non-negativity of all components in the solution, implies that $S(t), \, \, \, I(t)$ and $R(t)$ are all bounded. By the theory on the continuation of solutions for functional differential equations (see e.g., ^[26]), the boundedness of the solution also implies that the $T_m = \infty$, i.e., the solution exists globally.

Combining the above, we have actually established the well-posedness of the model (1.3)- (1.4) as stated in the following Theorem.

Theorem 2.1. For every $\phi = \left(\phi_{1}, \phi_{2}, \phi_{3}\right)^{T} \in C_+$, the model system (1.3)-(1.4) has a unique solution $(S(t), I(t), R(t))$ satisfying (2.1), which exists globally in $(0, \infty)$, is nonnegative and remains bounded. Moreover, the region

is both positively invariant and attractive for models (1.3)-(1.4).

3.   Disease free equilibrium and basic reproduction number

It is easy to see that the model (1.3) has the disease free equilibrium $E_0 = (S_0, 0, 0),$ where $S_0 = \Lambda/d$. To explore the stability of $E_0$ and identify the basic reproduction number of the model system (1.3)-(1.4), we linearize the model at $E_0$ to obtain

This is a linear ODE system (although the model (1.3)-(1.4) is a DDE system) with the coefficient matrix

which has two negative eigenvalues are $-d < 0, \, \, -(d+\alpha) < 0$, with the third eigenvalue $\beta P(0) S_0-(d+r+\epsilon), -(d+\alpha)$ being negative if and only if $\beta P(0) S_0 < (d+r+\epsilon)$. Thus $E_0$ is asymptotically stable if $\beta P(0)S_0 < (d+r+\epsilon)$ and unstable if $\beta P(0) S_0 > (d+r+\epsilon)$.

Actually, we can prove that $E_0$ is globally asymptotically stable when $\beta P(0)S_0 < (d+r+\epsilon)$. Indeed, when proving the boundedness of solutions, we have seen that

Choose a $\delta > 0$ sufficiently small such that $\beta P(0)(S_0 +\delta) < (d+r+\epsilon)$. Then, there is a $T > 0$ such that

By the conditions on $P(L)$, the second equation in (1.3) has the following linear ODE as an comparison equation from large $t$ from above

Note that all solutions converge to $0$ as $t \to \infty$ (since $\beta P(0)(S_0 +\delta) < (d+r+\epsilon)$). By a comparison argument, the $I(t)$ component of the solution to the model (1.3)-(1.4) also tends to $0$ as $t \to \infty$. This implies that the third equation in (1.3)-(1.4) has a limit equation $R'(t) = -(d+\alpha) R(t)$ which has the global convergence dynamics: $R(t) \to 0$ as $t \to \infty$; and the first equation in (1.3)-(1.4) has also a limit equation: $S'(t) = \Lambda -d S(t)$ which also have the global convergence dynamics: $S(t) \to \Lambda/d = S_0$ as $t \to \infty$. Now, by the theory of asymptotically autonomous systems (see, e.g., Castillo-Chavez et al. ^[27]), every solution with non-negative initial functions satisfies converges to $(S_0, 0, 0) = E_0$ as $t \to \infty$ when $\beta P(0)S_0 < (d+r+\epsilon)$.

The basic reproduction number of the model can be obtained from (3.1) by using the next generation approach. However, here it is more straightforward by tracking the duration of infection which is $1/(d+r+\epsilon)$, and transmission rate near the disease free equilibrium which is $\beta P(0)S_0$. Such a tracking immediately leads to the following formula for the basic reproduction number $R_0$:

Note that $\beta P(0)S_0 < (d+r+\epsilon)$ is equivalent to $R_0 < 1$.

When $R_0 > 1$ ($\beta P(0) S_0 > (d+r+\epsilon)$), $E_0$ becomes unstable and one would expect an endemic equilibrium. Note that because of the normality condition for the weights $k_j, \, \, j = 0, 1, \cdots, n$ in (1.2), possible endemic equilibrium $E^* = (S^{*}, I^{*}, R^*)$ is determined by the following equations

The first two equations in (3.4) can be rewritten as

$F(I)$ is decreasing since $\frac{r\alpha}{d+\alpha}-(d+r+\epsilon) < -(d+\epsilon) < 0$; while $G(I)$ is increasing because $P(L)$ is decreasing. Hence, (3.4) has a positive solution (unique) if and only if $G(0) < F(0)$, that is,

which is equivalent to $R_0 > 1$. Thus, we have proved that the model (1.3)-(1.4) has a unique endemic equilibrium $E^* = (S^{*}, I^{*}, R^*)$ if and only if $R_0 > 1$.

Summarizing the above analysis, we have actually proved the following Theorem.

Theorem 3.1. The disease-free equilibrium $E_0$ for the model (1.3)-(1.4) is globally asymptotically stable if $R_0 < 1$; and if $R_0 > 1$, $E_0$ becomes unstable and there occurs a unique endemic equilibrium $E^* = (S^*, I^*, R^*)$ determined by (3.4).

4.   Stability of the endemic equilibrium for a particular P(L(t))

In the preceding section, we have shown that when $R_0 > 1$, the model (1.3)-(1.4) has a unique endemic equilibrium $E^* = (S^{*}, I^{*}, R^*)$. This section is devoted to the stability of $E^*$. Unfortunately, for general $P(L(t))$, it does not seem to be possible to obtain any meaningful results on stability of $E^*$. Thus, in this section we consider the following particular severity measurement $L(t) = k_0I(t)+k_1I(t-\tau),$ which only looks at the weighed total infected cases at the present time $t$ and a previous time $t-\tau$, and the exponential decay function for $P(L)$: $P(L) = e^{-hL}$. The above particular choices result in the infection force $f(t) = \beta I(t) e^{-h(k_0I(t)+k_1I(t-\tau))}$ and accordingly reduce the model (1.3)-(1.4) to the following more concrete model system:

Now the basic reproduction number reduces to $R_0 = \frac{\beta \Lambda}{d(d+r+\epsilon)}$. Note that when $k_0 = 0$ and $\alpha = 0$, (4.1) further reduces to (1.5).

Assume $R_0 > 1$ so that the unique endemic equilibrium $E^* = \left(S^{*}, I^{*}, R^*\right)$ exists. The linearization of (4.1) at $E^*$ is given by

where $m = d+r+\epsilon$. From (4.2), we can derive the characteristic equation (CE) as

where

When $\tau = 0$, the transcendental equation (4.3) reduces to the following polynomial:

where

Straightforward calculation also shows that

By the Routh$-$Hurwitz criterion, all roots of the above cubic equation $F_1(\lambda, 0) = 0$ have negative real parts, implying that $E^* = (S^*, I^*, R^*)$ is locally asymptotically stable when $\tau = 0$.

Next, we explore the possibility for (4.3) to have a root with positive real part when $\tau$ is increased. This can occur only when a root of (4.3) crosses the purely imaginary axis in the complex plane from the left half to the right half when $\tau$ increases from zero. Noting that

for all $\tau \ge0$, the aforementioned crossing cannot happen at the origin of the complex plane. Thus, a crossing can only possibly occur in pairs of purely imaginary roots when $\tau$ is increased to some critical values.

To investigate this possibility, we plug $\lambda = i\omega$ (without loss of generality, we assume $\omega > 0$) into (4.3). Separating the real part and imaginary part of in the resulting equation $F_1(i\omega, \tau) = 0$ leads to

where

By eliminating the trigonometric functions, we obtain (4.5) the following equation:

Denoting $x = \omega^2$, we obtain

where

Thus, if $F_2(x) = 0$ has no positive root, then $E^*$ remains asymptotically stable for all $\tau > 0$. If $F_2(x) = 0$ has a positive root $x > 0$, then from (4.5), $\omega = \sqrt{x}$ would satisfy

which yields a sequence $\{\tau^n\}$, $n = 0, 1, 2, 3...$, of critical values for the delay parameter $\tau$, given by

To verify the transversality at these critical values, we first establish the following lemma and the detailed proof is shown in the Appendix.

Lemma 4.1. $\operatorname{sign}\left(\left.\operatorname{Re}\left(\frac{{\mathrm{d}} \lambda}{{\mathrm{d}} \tau}\right)\right|_{\tau = \tau^{n}}\right) = \operatorname{sign}\left(\left.\frac{{\mathrm{d}} F_2(x)}{{\mathrm{d}} x}\right|_{x = \omega^{2}}\right)$, where $\tau^{n}$ is given by (4.10).

We have seen that when $\tau = 0$, $E^*$ is asymptotically stable. Thus, the first critical value $\tau^0$ is of the most significance, as it is the smallest possible value for $\tau$ at which $E^*$ may lose its stability when $\tau$ passes it. In other words, we need to be mainly concerned with the case when there is a pair of eigenvalues crossing from the left complex plane to the right complex plane when $\tau$ is increased to pass $\tau^0$. Lemma 4.1 shows that only when the positive root $x$ of $F_2(x)$ satisfied $F'_2(x) > 0$, can such a crossing occur due to the increase of $\tau > 0$. Noticing that $F_2(x)$ is a cubic polynomial, it can have at most two positive roots at which $F_2'(x) > 0$. The crossing direction at $\pm i \sqrt{x}$ when $x$ is the largest positive root of $F_2(x) = 0$ (as $\tau$ passes a critical value given in (4.10)) is always from left to right; at $\pm i \sqrt{x}$ when $x$ is the second largest positive root of $F_2(x) = 0$ (if any) is from right to left. Thus, we have the following cases.

(C1) If $F_2(x) = 0$ has no positive root, then all roots of (4.3) have negative real parts.

(C2) If there is only one positive root $x_1$, then the crossing at $\pm \sqrt{x_1}$ rightward when $\tau$ increases to pass $\tau^0$; in this case, all roots of (4.3) have negative real parts for $\tau < \tau^0$, and (4.3) has a pair of root with positive real parts when $\tau > \tau^0$. Here $\tau^0$ is defined by (4.10) with $\omega = \sqrt{x_1}$).

(C3) If there are two positive roots $x_1 > x_{-} > 0$ at which $F'_2(x_1) > 0$ and $F'_2(x_{-}) < 0$ respectively. The crossing direction at $\pm i\sqrt{x_1}$ when $\tau$ passes $\tau = \tau_1^{n, +}$ is rightward, but at $\pm i\sqrt{x_{-}}$ the crossing direction as $\tau$ passes $\tau = \tau^{n, -}$, $n = 0, 1, 2...$ is leftward. Unlike case (C2), in this case, it is possible for the branch of the root $\lambda (\tau)$ that enters the right half of the complex plane through $\pm i \sqrt{x_1}$ as $\tau$ passes $\tau_1^{0, +}$, to return to the left half of the complex plane when $\tau$ further increases to $\tau^{0, -}$ (defined by (4.10) with $\omega_- = \sqrt{x_-}$).

(C4) If there are three positive roots: $x_1 > x_{-} > x_2$ at which $F_2'(x_i) > 0$, $i = 1, 2$ and $F_2'(x_{-}) < 0$ respectively, then the crossings at $\pm i\sqrt{x_1}$ and $\pm i\sqrt{x_2}$ when $\tau$ increases to pass $\tau = \tau_1^{n, +}$ and $\tau = \tau_2^{n, +}$ (defined by (4.10) with $\omega_1 = \sqrt{x_1}$ and $\omega_2 = \sqrt{x_2}$ respectively) are both rightward, but the crossing at $\pm i\sqrt{x_{-}}$ is leftward when $\tau$ increases to pass $\tau = \tau^{n, -}$ (defined by (4.10) with $\omega_- = \sqrt{x_-}$).

Combining the above summarization with the Hopf Bifurcation Theorem for delay differential equations, we accordingly proved the following theorem.

Theorem 4.1. Assuming $R_0 > 1$ so that $E^*$ exists. Then there can be four cases.

(i) If $F_2(x) = 0$ has no positive root, then $E^*$ is locally asymptotically stable for all $\tau\geq 0$.

(ii) If $F_2(x) = 0$ has only one positive root $x_{1}$ at which $F_2^\prime(x_1) > 0$, then $E^*$ is locally asymptotically stable for $\tau \in[0, \tau^{0})$ and unstable for $\tau > \tau^{0}$, where $\tau^{0}$ is given in (4.10) with $\omega = \sqrt{x_1}$. Moreover, system (4.1) undergoes Hopf bifurcation around $E^*$ at $\tau = \tau^0$.

(iii) If $F_2(x) = 0$ has two positive roots $x_{1}$ and $x_{-}$ at which $F_2^\prime(x_1) > 0$ and $F_2^\prime(x_{-}) < 0$ respectively, then $E^*$ is locally asymptotically stable for $\tau \in[0, \tau_1^{0, +})$ and becomes unstable when $\tau$ passes through $\tau = \tau_1^{0, +}$, where $\tau_1^{0, +}$ is given in (4.10) with $\omega = \sqrt{x_1}$. Moreover, system (4.1) undergoes Hopf bifurcation around $E^*$ at $\tau = \tau_1^{0, +}$.

(iv) If $F_2(x) = 0$ has three positive roots: $x_1 > x_{-} > x_2$ at which $F_2'(x_i) > 0$, $i = 1, 2$ and $F_2'(x_{-}) < 0$ respectively, then $E^*$ is asymptotically stable for $\tau \in\left[0, \tau^{0}\right)$ and becomes unstable when $\tau$ passes through $\tau = \tau^{0}$, where $\tau^{0} = \min \left\{\tau_1^{0, +}, \tau_2^{0, +}\right\}$ with

and

Moreover, system (4.1) undergoes Hopf bifurcation at $\tau = \tau^{0, +}_{ 1}$ and $\tau = \tau^{0, +}_2$.

We point out that for cases (ii)–(iv) in the above theorem, $E^*$ loses its stability when $\tau$ increases to $\tau^0$ or $\tau_1^{0, +}$; however addition, for case (iii) and (iv), $E^*$ may regain its stability at some $\tau > \tau^{0, -}$. We will numerically demonstrate this later.

Theorem 4.1 is a general conclusion, stated in terms of the delay $\tau$. To better understand the role of allocation of the weights $k_0$ and $k_1$ in detail, we can do further specific analysis on positive roots of the cubic equation $F_2(x) = 0$ below, in terms of $k_0$ and $k_1$. To this end, we denote

where $q_j(k_1), \, \, j = 0, 1, 2,$ are give in (4.8). It can be readily seen that $q_0(k_1)$, $q_1(k_1)$ and $q_2(k_1)$ are decreasing in $k_1$ and each of them has a unique positive root $\tilde{k}_{1i}, \, \, \, i = 0, 1, 2$. From (4.8), plugging the formulas for $r_j, \, \, \, j = 0, 1, 2$, and after some tedious calculations (omitted here), these roots are given by

with $q_1(0), \, \, q_1(1), \, \, q_2(0)$ and $q_2(1)$ explicitly given by

It can be readily seen that $q_0(k_1)$, $q_1(k_1)$ and $q_2(k_1)$ are decreasing in $k_1$. In addition, $q_i(\tilde{k}_{1i}) = 0$ ($i = 0, 1, 2$). Here tedious calculations for the above formulas are omitted.

By the definition of $\tilde{k}_{1i}, \, \, i = 0, 1, 2$, we immediately have the following.

Lemma 4.2. $q_i(k_1) > 0$ iff $k_1 < \tilde{k}_{1i}$, where $i = 0, 1, 2$.

According to reference [17] (Lemma A.2 and A.3) or references [28,29], as well as Lemma 2 above, we obtain the following results.

Corollary 4.1. Assume that $R_0 > 1$.

(I) Case (ii) in Theorem 4.1 occurs if and only if one of the following conditions is satisfied

1) $\Delta(k_1)\geq 0$ and $\tilde{k}_{10} < k_1$;

2) $\Delta(k_1) < 0, \, \, \tilde{k}_{10} < k_1 < \min\{\tilde{k}_{11}, \tilde{k}_{12}\}$;

3) $\Delta(k_1) < 0, \, \, \max\{\tilde{k}_{11}, \tilde{k}_{10}\} < k_1$;

4) $\tilde{k}_{10} < k_1 = \tilde{k}_{11}$ or $\tilde{k}_{11} < k_1 = \tilde{k}_{10}$;

(II) Case (iii) in Theorem 4.1 occurs if and only if one of the following conditions is satisfied

1) $\Delta(k_1) < 0, \, \, \tilde{k}_{12} < k_1 < \min\{\tilde{k}_{10}, \tilde{k}_{11}\}$;

2) $\Delta(k_1) < 0, \, \, \tilde{k}_{11} < k_1 < \tilde{k}_{10}$;

(III) Case (iv) occurs if and only if there holds (v): $\Delta(k_1) < 0$, $\max\{\tilde{k}_{10}, \tilde{k}_{12}\} < k_1 < \tilde{k}_{11}$;

(IV) Case (i) occurs if and only if none of (I)–(III) is satisfied. Moreover, Case (i) holds if $k_1\leq\min\{\tilde{k}_{10}, \tilde{k}_{11}, \tilde{k}_{12}\}$ (i.e. $q_i(k_1)\geq0, i = 0, 1, 2$).

Particularly, if $k_1 > \max\{\tilde{k}_{10}, \tilde{k}_{11}\}$ holds, then one of cases (I)-1) and (I)-3) must occur. Based on this observation, we obtain

Corollary 4.2. If $k_1 > \max\{\tilde{k}_{10}, \tilde{k}_{11}\}$ (i.e. $q_i(k_1) < 0, i = 0, 1$), then Case (ii) in Theorem 4.1 occurs.

From Corollaries 4.1 and 4.2, a necessary condition for Case (ii) of Theorem 4.1 is $\tilde{k}_{10}\leq k_1$, and a necessary condition for Case (iii) of Theorem 4.1 is $\tilde{k}_{10} > k_1$. Interestingly, according to (IV) in Corollary 4.1, there will be no bifurcation at all provided that the weight $k_1$ is small ($k_1\leq\min\{\tilde{k}_{10}, \tilde{k}_{11}, \tilde{k}_{12}\}$). In other words, restricting the weight to the previous surveillance $I(t-\tau)$ to be sufficiently small can avoid Hopf bifurcation (sustained oscillation). From Corollary 4.2, on the other hand, a relatively large weight $k_1$ will allow $\tau$ to destroy the stability of $E^*$ at some critical value. In Section 6, we will numerically illustrate the impact of $k_1 = 1-k_0$ in terms of these two corollaries.

5.   Special case $\alpha=0$, $\epsilon=0$: a comparison to related work

In order to compare our results with those in paper ^[16], let us also consider the infection force $f(t) = \beta I(t) e^{-h(k_0I(t)+k_1I(t-\tau))}$ and let $\alpha = 0$, $\epsilon = 0$ in the model (4.1), reducing the model (4.1) to the following model system

The case $k_0 = 0$ (hence $k_1 = 1$) leads to the model (1) in reference [16] (i.e., Eq.(1.5) in the introduction of this paper), which has been thoroughly analyzed. Note that $k_0 = 0$ means that only infected cases in the previous time $t-\tau$ are used to measure the severity $L$. Here in this section, we allow $k_0\in (0, 1]$ (hence $k_1 \in [0, 1)$ and hence, the severity $L$ is measured by the weighted total number of infected cases at the present time $t$ and the past time $t-\tau$. This seems to be a more realistic and reasonable choice as it allows public health agents to put different weight on the disease information at different time; in the mean time, this also constitutes a good demonstration of flexibility of our model framework in accommodating a wide range of scenarios. By analyzing (5.1) in more details, we can reveal impact of the weights $k_0$ and $k_1$, together with the time delay $\tau$ (see a discussion in Section 7).

The equilibria structure and the basic reproduction number $R_0$ for (5.1) remain the same as in reference [16] for (1.5), due to the complementary constraint $k_0+k_1 = 1$. The first two equations of (5.1) form a decoupled sub-system:

Let $m = d+r$. The CE at $E^*$ for (5.2) becomes

in which

with

To be consistent with Section 4, we use the notations here:

and

Lemma 4.1 also holds for this model, but now function $F_2(x)$ becomes a quadratic function (instead of a cubic function):

where

Let

When $q_0(k_1) < 0$, $F_2(x)$ has a unique positive root: $x_1 = -\frac {q_2(k_1)}{2}+\frac {\sqrt {\Delta(k_1)}}{2}$ at which $F_2'(x_1) > 0$. Thus, $E^*$ loses its stability for all $\tau > \tau^0$ to a periodic solution through Hopf bifurcation.

When $q_0(k_1) > 0, \, q_2(k_1) < 0$ and $\Delta(k_1) > 0$, there are two positive roots $x_1 = -\frac {q_2(k_1)}{2}+\frac {\sqrt {\Delta(k_1)}}{2}$ and $x_{-} = -\frac {q_2(k_1)}{2}-\frac {\sqrt {\Delta(k_1)}}{2}$ such that $F_2'(x_1) > 0$ and $F_2'(x_{-}) < 0$. Then (5.3) has two pairs of simple purely imaginary roots $\pm i\omega_{+}$, $\pm i\omega_{-}$ at $\tau^{n, +}$, $\tau^{n, -}$, $n = 0, 1, 2...$, where $\omega_{+} = \sqrt{x_1}$, $\omega_{-} = \sqrt{x_{-}}$ and $\tau^{n, +}$, $\tau^{n, -}$ are defined by the same form in (4.10) with $\omega$ given by $\omega_{+} = \sqrt{x_1}$, $\omega_{-} = \sqrt{x_{-}}$ respectively. According to Lemma 4.1, we have

Hence, at each $\tau = \tau^{n, +}$, a pair of roots crosses the imaginary axis at $\pm i \omega_{+}$ into the right half-plane, and at each $\tau = \tau^{n, -}$, a pair of roots crosses at $\pm i \omega_{-}$ back into the left half-plane. These observations together with the fact that all roots of (5.3) are on the left half of the complex plane when $\tau = 0$, imply that $\tau^{0, +} < \tau^{0, -}$. However, the orders of other numbers in the two sequences $\tau = \tau^{n, +}$ and $\tau = \tau^{n, -}$ do not seem to have a pattern and seem to be complicated.

For example, for the first couple of critical values in the two sequences, the following two cases are possible:

(A) $\tau^{0, +} < \tau^{0, -} < \tau^{1, +}$. In this case, all roots of (5.3) have negative real parts when $\tau \in [0, \tau^{0, +})$; one pair of roots of (5.3) have positive real parts when $\tau \in (\tau^{0, +}, \tau^{0, -})$ and that pair will return to the left complex plane again when $\tau \in (\tau^{0, -}, \tau^{1, +})$.

(B) $\tau^{0, +} < \tau^{1, +} < \tau^{0, -}$. In this case, all roots of (5.3) have negative real parts when $\tau \in [0, \tau^{0, +})$; one pair of roots of (5.3) have positive real parts when $\tau \in (\tau^{0, +}, \tau^{1, +})$ and there will be two pairs of roots on the right complex plane when $\tau \in (\tau^{1, +}, \tau^{0, -})$.

Summarizing the above, we have the following:

Theorem 5.1. Assuming $R_0 > 1$, then

(i) If $q_0(k_1) < 0$, or $q_0(k_1) = 0$ but $q_2(k_1) < 0$, the endemic equilibrium $E^*$ of system (5.1) is locally asymptotically stable for $\tau\in[0, \tau^0)$ and unstable for all $\tau > \tau^0$; moreover, system (5.1) undergoes Hopf bifurcation when $\tau$ passes $\tau = \tau^{n}, n = 0, 1, 2 \ldots,$ where $\tau^n, \, \, \, n = 0, 1, 2, \cdots$ are defined by (4.10) with $\omega = \sqrt{x_{1}}$.

(ii) If $q_0(k_1) > 0$, $q_2(k_1) < 0$ and $\Delta(k_1) > 0$, then there are two sequences$\{\tau^{n, +}\}$ and $\{\tau^{n, -}\}$, still defined by ((4.10)) but with $\omega$ being$\sqrt{x_1}$ and $\sqrt{x_-}$ respectively, such that$E^*$ is locally asymptotically stable when $\tau\in[0, \tau^{0, +})$, and it loses its stability when $\tau$ increases to pass $\tau^{0, +}$ through Hopf bifurcation. It is possible for $E^*$ to regain its stability for some $\tau > \tau^{0, -}$. Moreover, system (5.1) undergoes Hopf bifurcation at each $\bar{\tau}\in \{\tau^{n, +}\}\cup \{\tau^{n, -}\}$.

(iii) For other cases of $q_0(k_1)$ and $q_2(k_1)$, $E^*$ is asymptotically stable for all $\tau\ge 0$.

To obtain more details on how the weight parameter $k_1$ affects the conditions in Theorem 5.1, we need to explore the signs of $q_0(k_1)$ and $q_2(k_1)$ given by (5.5). Both $q_0(k_1)$ and $q_2(k_1)$ are decreasing in $k_1$ and each has a unique root $\tilde{k}_{1i} > 0$ which can be expressed in terms of the model parameters by

where

Thus, a version of Lemma 4.2 also holds here: $q_{i}(k_1) > 0$ iff $k_1 < \tilde{k}_{1i}, \, \, \, i = 0, 2$. Moreover, if $\tilde{k}_{12}\leq\tilde{k}_{10}$, then there exists

such that $\Delta(k_{ic}) = 0$ ($i = 1, 2$), and $\Delta(k_1) > 0$ when $k_1 > k_{1c}$ or $k_1 < k_{2c}$.

By symmetry, using the same arguments for Theorem 4.1 and Corollary 4.1, we have the following more specific results for Theorem 5.1, in terms of $k_1$:

Since Lemma 4.2 also holds, then we get more specific results as follows:

Corollary 5.1. Assuming $R_0 > 1$, then

(i) If $k_1 > \tilde{k}_{10}$ or $\tilde{k}_{12} < k_1 = \tilde{k}_{10}$, the conclusion of Theorem 5.1-(i) holds.

(ii) If $\tilde{k}_{12} < k_1 < \tilde{k}_{10}$ and $\Delta(k_1) > 0$ (i.e., $\tilde{k}_{12} < k_{1c} < k_1 < \tilde{k}_{10}$), then the conclusion of Theorem 5.1-(ii) holds.

(iii) For other range of $k_1$, particularly when $k_1\leq \min\{\tilde{k}_{12}, \tilde{k}_{10}\}$, the conclusion of Theorem 5.1-(iii) holds.

As the above results show, when $k_1$ is small to a certain extent (i.e., $k_1\leq \min\{\tilde{k}_{12}, \tilde{k}_{10}\}$), the stability of $E^*$ can be ensured whatever $\tau > 0$. When $k_1$ is large to a certain extent (i.e., $k_1 > \tilde{k}_{10}$), the system (5.1) will undergo Hopf bifurcations at a sequence $\{\tau^{n}\}$ as $\tau$ increases. We point out that when $k_1 = 1$ ($k_0 = 0$), the condition $1 = k_1 > \tilde{k}_{10}$ in (i) of Corollary 5.1 reduces to $\frac{m}{2S^*dh} < 1$ which is equivalent to $\frac{dh}{\beta}\exp(\frac{\Lambda h}{m}-1) < 1$ by the properties of the Lambert W function ^[30], and thus, is consistent with the conclusion of Case1: (H1) in reference [16]. Corollary 5.1-(ii) corresponds to the conclusion of Case2: (H2) in reference [16]. Because of the existence of two sequences of critical values $\{\tau^{n, +}\}$ and $\{\tau^{n, -}\}$ which may allow different order patterns (see, e.g., (A) and (B)), the above analysis indicates that the stability switches as $\tau$ increase for this case can be more complicated than stated in Proposition 4 in reference [16]; for example, there may be multiple switches between the stability and instability of $E^*$ as $\tau$ increases, and this will be demonstrated numerically in the next section.

6.   Numeric exploration of multiple switches of stability

From the results in Sections 4 and 5, we have seen that within certain ranges of the parameter $k_1$, there may be one, two or even three sequences of critical values for the delay parameter $\tau$ at which Hopf bifurcations occur. Among these critical values (if any), the smallest $\tau^{0, +}$ is most important because the endemic equilibrium $E^*$ is deemed to lose its stability when $\tau$ increases to it. However, the effect of the remaining critical values (if any) depends on the order they appear. Particularly when the $F_2$ function in Section 4 or 5 has more than one positive root, it is generally very difficult (if not impossible) to determine whether there is a pattern for the order of those critical values in the two or three sequences of critical values for $\tau$ related to the two or three positive roots of $F_2$. Below we provide some numeric examples to show that such orders can be different depending on the values of the model parameter. Our numerical diagrams are constructed mainly by using the DDE Biftool package ^[31,32,33].

We begin with (5.1). For convenience of comparison with the results in reference [16], we choose the same values of the parameters for (5.1) as used in reference [16] as below:

By numerical computations, we obtain $E^*\approx(0.62946, 0.24703, 0.12351)$ and $R_0 = 3.3$, together with the three threshold values of $k_1$: $\tilde{k}_{12}\approx0.8735925, \, \, \, k_{1c}\approx0.89673437$ and $k_{2c} < 0$. Then $k_1 > k_{1c}$ (resp. $k_1\in [0, k_{1c})$) implies $\Delta(k_1) > 0$ (resp. $\Delta(k_1) < 0$).

The related conclusions corresponding to Corollary 5.1 are illustrated in the $k_1-\tau$ plane in Figure 1 and are restated below:

a) When $k_1 < k_{1c}$, no Hopf bifurcation occurs when $\tau$ increases, and $E^*$ remains stable for all $\tau\ge0$ (Figure 1(a));

b) For $k_{1c} < k_1 < \tilde{k}_{10}$ (see the zoom-in figure in Figure 1(b)), there are two sequences of critical values beginning with $\tau^{0, +}$ and $\tau^{0, -}$, respectively. The stability of $E^*$ switches as $\tau$ increases from zero: $E^*$ first loses stability at $\tau^{0, +}$ and regains the stability at $\tau^{0, -}$;

c) For $\tilde{k}_{10} < k_1\leq 1$, there is only one sequence of critical values beginning with $\tau^{0, +}$; $E^*$ loses stability at $\tau = \tau^{0, +}$ and remains unstable for all $\tau > \tau^{0, +}$ (see the zoom-in figure Figure 1(b));

d) When $\tau < 11.7575$, $E^*$ remains stable regardless of $k_1 \in [0, 1]$ (see Figure 1(a)).

Especially, when $k_1 = 0.8968\in(k_{1c}, \tilde{k}_{10})$, we can calculate the first few critical values to obtain the following alternating order:

with $\omega_{+} = 0.0597$, $\omega_{-} = 0.0326$. Accordingly, when $\tau$ increases from $\tau = 0$ to pass $\tau^{2, +}$, the stability/instability $E^*$ switches at least five times: losing stability at $\tau^{i, +}, \, \, i = 0, 1, 2$, and regain stability at $\tau^{i, -}, \, \, i = 0, 1$. See Figure 2.

When $k_1 = 1$, bifurcation also occurs at $\tau^n$ with $\tau^0\approx11.7575$ and no re-stabilization of $E^*$ occurs when $\tau > \tau^0$, which is in line with the conclusion and Figure 1 in reference [16].

If we fix $\tau$ at a given value, we can also use our results to numerically demonstrate the impact of $k_1$. This corresponds to a scenario that current data and the data of fixed time ago, are available and we want to explore how the weight allocation weight $(k_0, k_1)$ will affect the disease dynamics.

Next, we look at the more general (4.1). Fix the other parameters in (4.1) as below

Then by numeric computation we obtain $E^*\approx(0.65043, 0.20657, 0.03972)$, $R_0 = 2.86$, $\tilde{k}_{10}\approx0.93364$, $\tilde{k}_{11}\approx0.88864$, $\tilde{k}_{12}\approx1.27804$ and $\Delta(k_c) = 0$ with $k_c = 0.929195195$. Note that $k_1 > k_c$ (resp. $k_1\in[0, k_c)$) implies $\Delta(k_1) < 0$ (resp. $\Delta(k_1) > 0$). By Corollary 4.1, we have the following conclusions for (4.1) which are similar to a)–d) above for (5.1):

(a When $0 < k_1 < k_c$, no Hopf bifurcation occurs when $\tau$ increases, and $E^*$ remains stable for all $\tau\ge 0$;

(b When $k_c < k_1 < \tilde{k}_{10}$, there are two sequences of critical values for $\tau$ starting with $\tau^{0, +}$ and $\tau^{0, -}$, respectively. There are multiple switches of the stability/instability of $E^*$ as $\tau$ increases from zero: $E^*$ first loses stability at $\tau^{0, +}$ and regains the stability at $\tau^{0, -}$; and then loses stability again at $\tau^{1, +}$;

(c When $\tilde{k}_{10} < k_1\leq 1$, there is only one sequence of critical values for $\tau$ starting with $\tau^{0, +}$, and $E^*$ loses stability at $\tau = \tau^{0, +}$ remains unstable for all $\tau > \tau^{0, +}$;

(d If $\tau < 12.45 = \tau^{0, +}(1)$, then $E^*$ remains stable for all $k_1 \in [0, 1]$.

The first bifurcation branch in terms of $k_1$ (i.e., $\tau^{0, +}$ and $\tau^{0, -}$) is demonstrated in Figure 3, which visually reflects the above conclusions (a–(d.

Although the first bifurcation branches in Figure 3 and the zoom-in Figure 1(a) are quantitatively the same, the order of subsequential branches can be different. To see this, we choose $k_1 = 0.931\in(k_c, \tilde{k}_{10})$. Then there are two positive roots of $F_2(x) = 0$: $x_1 = 0.0099$ with $F_2'(x_1) > 0$ and $x_{-} = 0.0021$ with $F_2'(x_{-}) < 0$; and by numerical calculations, we obtain $\omega_{1} = 0.0997$, $\omega_{-} = 0.0453$, and accordingly

which differs from the alternating order in (6.2). The switches of stability/instability of $E^*$ is shown in Figure 4(a), which has different pattern of zeros from that in Figure 2(a). For the above parameter values, Figure 6(a), (b) illustrate the solutions of (4.1) with initial value $(0.7, 0.2, 0.04)$: for $\tau = 15 < \tau_1^{0, +}$, the solution tends to $E^*$; when $\tau = 30\in (\tau_1^{0, +}, \tau^{0, -})$, the solution converges to a periodic orbit surrounding $E^*$.

If fixing $k_1 = 0.94\in[\tilde{k}_{10}, 1]$, then $F_2(x) = 0$ has a unique positive root $x_1 = 0.0163 > 0$ at which $F_2'(x_1) > 0$. For this case, there is only one sequence of critical values $\{ \tau^{n} \}$ with $\tau^{0} = 20.82$, $\tau^{1} = 69.97$, $\tau^{2} = 119.13$ and $\tau^{3} = 168.28$, shown in Figure 4(b). Noticing that $0.94$ and $0.931$ are very close, this indicates that the occurrence of Hopf bifurcation from the endemic equilibrium $E^*$ with respect to the two parameters and the path of bifurcation can be sensitive and complicated. This phenomenon was not explored in reference [16] since the model in reference [16] has $k_1 = 1$.

7.   Conclusions and discussion

In this paper, we have provided a framework from a new perspective to look at infection forces in modelling infectious disease transmission dynamics. This new angle was motivated by the fear effect in predator-prey interactions as well as by the effects of all those non-pharmaceutical control measures implemented around the world during the ongoing Covid-19 pandemics. Such a framework not only can accommodate various infection force functions in existing infectious disease models, but can also suggest new forms of infection forces, particularly infection forces that depend on the past disease surveillances. We have demonstrated this new framework by incorporating, into a classic SIRS model, an infection force that depends on both the disease surveillance at current time $t$ and the disease surveillance at a past time $t-\tau$, with each being given a weight ($k_0$ and $k_1$, respectively). The resulting new model (4.1) is a system of delay differential equations, which includes some existing models as special cases. Here the two new parameters $\tau$ and $k_1$ ($k_0 = 1-k_1$) can be considered as kind of strategy parameters, reflecting the impact of the current and past disease surveillances on disease control measures and human behaviours.

We have performed a standard analysis on the new model (4.1), including confirming well-posedness of the model, identifying the basic reproduction number $R_0$ for the model, establishing the threshold dynamics for the model in terms of $R_0$, and exploring the stability of and bifurcation from $E^*$ with respect to the two key parameters $\tau$ and $k_1$. Our results show that although the basic reproduction number of the new model (featured by $\tau$ and $k_1$ compared to existing models) does not affect $R_0$ and the threshold between disease-free and endemic, they do affect the disease dynamics (long-time disease patterns) when it becomes endemic. More specifically, we have found that there are ranges on $\tau$ and $k_1$ for which the model demonstrates convergent (to $E^*$) endemic pattern and oscillatory endemic pattern (due to Hopf bifurcations) respectively. Most interestingly, we have observed that the path of Hopf bifurcations when $\tau$ increases can be different when other model parameters are at different values, as illustrated by the first few critical values of $\tau$ in (6.2) and (6.4). Note that the model in the recent work ^[16] is a special case of our reduced model (4.1) where $k_1$ is fixed at 1 (hence $k_0 = 0$). Hence, even our reduced model (4.1) offers more flexibility in the use of disease surveillance data: depending on how recent the datum is (e.g., the value of $\tau$), its impact on the long term disease dynamics (through $k_1$) could be different. The results in Section 5, particularly the numerical results provided in Section 6 clearly show the role $k_1 \in [0, 1]$ can play and the difference it can make. These results on the impact of $k_1$ and $\tau$ on the long time dynamics of our models imply that the long-term effect of control measures changes when the delay in reporting of infection increases; therefore, they actually have demonstrated that timely reporting is important to management of a disease.

We point out that our results on Hopf bifurcations are only local. It is possible to carry out a global bifurcation analysis in the case that the equation $F_2(x) = 0$ only has one positive root, generating only one sequence of critical values for $\tau$. We choose not to do such a global bifurcation in this already lengthy paper since the main purpose is to demonstrate our new framework of modelling infection forces in transmission dynamics when non-pharmaceutical effects are considered. When $F_2(x) = 0$ has more than one positive root, the dynamics are even richer since multiple switches of stability/instability of endemic equilibrium may happen.

In reality, unavoidably there is a delay in the availability of disease surveillance at any given time. This is because it takes time to diagnose/test and report cases and collect the data during an epidemic. This fact forces one to use the past disease surveillances in measuring the severity of the epidemic, and hence justifies the form (1.1) and (1.2) for the severity function $L(t)$. Moreover, in addition to the number of infected cases $I(t)$, the rate of change of $I(t)$ (i.e., $I'(t)$) will also have an impact on the practically susceptible population (i.e., $S_p = P S$) through the proportion function $P(L)$: positive $I'(t)$ would increase $L$ (decrease $P(L)$) and negative $I'(t)$ would decrease $L$ (include $P(L)$). This idea was explored in Xiao et al. ^[34] where the authors proposed a media-impact function of the form $\beta_0e^{M(t)}$, where

Incorporating such a dependence on $I'(t)$ would significantly make the model more challenging to analyze, as seen in reference [34]. A natural alternative way to measure the rate of change is to look at an approximation of $I'(t)$ by the difference $\Delta I(t) = I(t)-I(t-\sigma)$ which has recently been used in Li et al. ^[35]. Note that the dependence of $L(t)$ on $\Delta I(t)$ can be absorbed into the form of (1.2). Therefore, the form (1.2) and the more general for (1.1) deserve more investigations in association with various proportion functions $P(L)$ (e.g., those in (A) and (B)) for infection disease models of various types.

Acknowledgments

The authors are grateful to the two anonymous referees for their valuable comments which have led to a substantial improvement in the presentation of the paper. This research is partially supported by NSERC of Canada.

In memory of Stephen Gourley

From the second author: I was so shocked by Stephen's passing away. I could not believe it and I refused to believe it for quite a while. Stephen was a great friend of mine for more than 20 years. As a friend, I enjoyed every meeting with him, every chat with him, every hiking with him, and every meal with him. As a collaborator, I shared many common research interests with him and I learnt a lot from him in our collaborations. Our friendship and collaborations later extended to some of my former students and postdocs - he was so nice and supportive to them, and I was very impressed by and grateful for his guidance and mentorship to these young people. Stephen was a true gentleman - friendly, caring, considerate, humorous, and stimulating; and hosting his visits was the easiest and most pleasant thing. Stephen will be remembered by me, and I am sure, by many others as well.

Conflict of interest

The authors declare that no competing interests exist.

Appendix

This appendix includes the detailed proof of the Lemma 1 stated in the present paper.

Proof. Take derivative of (4.3) about $\tau$:

then

This together with $F_1(i\omega, \tau^n) = 0$ further leads to

and

Therefore,