Mutual contacts among individuals exist in the real world. It is necessary to incorporate the interaction among infectious and susceptible individuals when studying the epidemic spread. Many approaches (based on nodes, pairwise, edges and so on) have been established and used to investigate the interaction of individuals in the network [1–5]. Contact tracing is a useful strategy in epidemic surveillance [6, 7]. The close contacts of an infected individual can be identified by contact tracing and have the chance to be infected. Screening and isolating these individuals is an effective way of identifying potentially infected individuals. By tracing the contacts of infected individuals, the identified infected individuals can be isolated and treated, and the transmission chains of epidemic can be detected and stopped so as to mitigate the spread of the epidemic.

Many mathematical models are developed to investigate the effect of contact tracing in controlling the spread of infectious diseases from the individual-based models (IBMs) to the ordinary differential equations (ODEs) at the population level. Kretzschmar et al. (1996) used IBMs to study the spread of two sexually transmitted diseases through simulations and found that contact tracing for targeting the highly sexually active core group is very effective in reducing the prevalence [8]. Kasaie et al. (2014) investigated the impact of household contact tracing for tuberculosis in a moderate-burden setting [9]. Peak et al. (2017) used IBMs to compare the effectiveness of contact tracing for Ebola, influenza, and SARS in simulations [10]. To obtain more analytic results on contact tracing, mean-field and pairwise approximation models are developed [11–13]. Eames (2007) utilized pairwise models on networks to compare the effects of recursive tracing and one-step tracing [14]. Heijne et al. (2010) applied pair approximation models to predict the impact of contact tracing on chlamydia prevalence [15]. Barlow (2020) studied the effect of contact tracing from the view of branching processes and phenomenological approaches [16].

In addition, the awareness of infection among infected individuals plays a key role in the spread of slowly progressive diseases. Many researchers focused on studying the roles of the public awareness among susceptible individuals in the spreading of diseases, and showed that raising awareness of protection among susceptible individuals can effectively contain epidemic spread [17–21]. However, to examine the effects of the infection awareness among infected individuals on the epidemic theoretically is still a great challenge. Some previous studies indicated that asymptomatic or mild infections are likely to cause outbreaks and rapid transmission in communities [22–24]. Even rapid contact tracing and antigen testing are effective in control the spread of infectious diseases and can be widely used, large numbers of infected individuals are unaware of their infection so that the epidemic spread is uncertain [25–27]. Considering the complexity of awareness process, we define the awareness process of infectious status, which starts from being infected to the time when realizing being infected, as a non-Markovian process in this study. The reason is that recent empirical observations show that non-Markovian transmission process, which refers to the spreading processes with non-exponential distributions, is more realistic [28–30]. Kiss et al. (2015) presented pairwise network models with non-Markovian recovery period to provide links between non-Markovian dynamics and pairwise delay differential systems [31]. Sherborne et al. analyzed the characteristics of an edge-based model with arbitrary transmission and recovery process [32]. Li et al. (2018) considered the fixed infectious period in an SIR pairwise model [33]. The existence of non-exponential distribution shows richer dynamics of epidemic spread [34–38], where the epidemic curves are different when the infectious periods have different distributions even with the same mean [39, 40].

In this work, we incorporate the non-Markovian process into pairwise models to consider the impacts of contact tracing and the awareness of infection status on the epidemic spread. Our work has following contributions. First, we derive the basic reproduction number and the contact tracing threshold for the non-Markovian pairwise model, and prove that the basic reproduction number will be less than one when the contact tracing rate is larger than the threshold. Second, we investigate the impacts of contact tracing and the awareness of infectious status on epidemic spreading. The work is organized as follows. We develop pairwise models incorporating contact tracing and non-Markovian awareness process for epidemic networks in Section 2. In Section 3, we analyze the positive invariance of the model, derive the basic reproduction number, and obtain the final epidemic size relation by rigorous analysis. Section 4 provides the theoretical results on the threshold of effective tracing rate and the effects of non-Markovian awareness processes. Numerical simulations are given in Section 5. Discussion is presented in Section 6.

We develop a pairwise network model with continuous non-Markovian awareness process to investigate the impacts of contact tracing on epidemics. Five classes are included in this model, namely, susceptible individuals (S), infected individuals (I), traced infected individuals (Q _c ), infected individuals who are aware of their own infectious status and treated (Q _a ), and recovered individuals (R). The susceptible individuals become infectious when they are infected. The infected individuals can be divided into three groups, that is, those who will be traced at the rate, τ, those who recover naturally at the rate, α, and those who are aware of infection and go to the hospital initiatively in a non-Markovian process. We consider the transition process from the state I to the state Q _a as the awareness process of infection among infected individuals, which can be described by a random variable, X. The probability density function of X is f(s). The corresponding cumulative function is F(s) and the survival function is ξ(s) = 1 − F(s). ϕ(s) represents the density of nodes from state I to state Q _a with respect to the age of s at the initial time. Furthermore, all individuals in classes Q _c and Q _a will be quarantined and treated. The infected individuals who are successfully treated will enter the recovered class at the rate, γ. The total number of individuals is N, i.e., [S] + [I] + [Q _c ] + [Q _a ] + [R] = N, where [A] is the number of individuals in state A. We use the notation [AB] to represent the expected number of links that connect a node in state A and a node in state B. Let m be the average degree in this network model.

The schematic diagram is depicted in Figure 1. Let i(t, a) be the density of the infected nodes with infection age a at time t, then the total number of infected nodes at time t is [ I ] ( t ) = ∫ 0 ∞ i ( t , a ) d a . Due to the McKendrick-von Foerster equation [41], the changes of i(t, a) can be given by the following partial differential equation ∂ i t , a ∂ t + ∂ i t , a ∂ a = − τ + α + h a i t , a , (1) in which h(a) = −ξ′(a)/ξ(a) = f(a)/ξ(a), i(t, 0) = β[SI](t), and i(0, a) = ϕ(a). Here, we assume that lim _a→∞ ϕ(a) = 0 according to the biological feasibility. From Eq. 1 and the boundary conditions of i(t, a), we deduce the representation of i(t, a) according to the work in [37] as follows. i t , a = β S I t − a e − ∫ t − a t τ + α d s e − ∫ 0 a h b d b , i f t > a . ϕ a − t e − ∫ 0 t τ + α d s e − ∫ a − t a h b d b , i f t ≤ a . (2) Then, using the properties of h(a), we obtain that d d t I t = ∫ 0 ∞ ∂ ∂ t i t , b d b = ∫ 0 ∞ − τ i t , b − α i t , b − h b i t , b − ∂ ∂ t i t , b d b = β S I t − τ + α I t − ∫ 0 ∞ h b i t , b d b = β S I t − τ + α I t − ∫ 0 t β S I t − s f s e τ + α s d s − ∫ t ∞ ϕ s − t f s e τ + α t ξ s − t d s . (3) Similarly, the expected number of links between S and I at time t is [ S I ] ( t ) = ∫ 0 ∞ S i ( t , a ) d a . Considering the removal of S-I links, we describe the deviation of Si(t, a) by the following partial differential equation. ∂ S i t , a ∂ t + ∂ S i t , a ∂ a = − β I S i t , a − β + τ + α + h a S i t , a , (4) where Si(t, 0) = β[SSI](t) and S i ( 0 , a ) ≈ m N [ S ] ( 0 ) ϕ ( a ) . To break higher order moments, we use the closure approximation formula [ A X B ] ( t ) = m − 1 m [ A X ] ( t ) [ X B ] ( t ) [ X ] ( t ) in [42]. Because d d t [ S I ] ( t ) = ∫ 0 ∞ S i ( t , a ) d a , we find d d t S I t = − ∫ 0 ∞ β m − 1 m S I t S t + β + τ + α + h a S i t , a d a − ∫ 0 ∞ ∂ ∂ a S i t , a d a = β 1 − m m S I 2 t S t − β + τ + α S I t − ∫ 0 ∞ h a S i t , a d a + β m − 1 m S S t S I t S t . (5) Solving the Eq. 4, we obtain the representation of Si(t, a) according to the work in [37]. S i t , a = β m − 1 m S S t − a S I t − a S t − a e − ∫ t − a t β m − 1 m S I s S s + β + τ + α d s e − ∫ 0 a h b d b , i f t > a ; m N S 0 ϕ a − t e − ∫ 0 t β m − 1 m S I s S s + β + τ + α d s e − ∫ a − t a h b d b , i f t ≤ a . (6) Applying the Eqs 5, 6, we have d d t S I t = β m − 1 m S S t S I t S t − β m − 1 m S I t S I t S t − β + τ + α S I t − ∫ 0 t β m − 1 S S t − s S I t − s f s m S t − s e ∫ t − s t β m − 1 m S I u S u + β + τ + α d u d s − ∫ t ∞ m S 0 ϕ s − t f s N e ∫ − 0 t β m − 1 m S I u S u + β + τ + α d u ξ s − t d s . (7) Therefore, the resulting pairwise model can be constructed with integro-differential terms in Model (8) as follows d d t S t = − β S I t , d d t I t = β S I t − τ + α I t − ∫ 0 t β S I t − s f s e τ + α s d s − ∫ t ∞ ϕ s − t f s e τ + α t ξ s − t d s , d d t S S t = − 2 β m − 1 m S S t S I t S t , d d t S I t = β m − 1 m S S t S I t S t − β m − 1 m S I t S I t S t − β + τ + α S I t − ∫ 0 t β m − 1 S S t − s S I t − s f s m S t − s e ∫ t − s t β m − 1 m S I u S u + β + τ + α d u d s − ∫ t ∞ m S 0 ϕ s − t f s N e ∫ 0 t β m − 1 m S I u S u + β + τ + α d u ξ s − t d s , d d t Q c t = τ I t − γ Q c t d d t Q a t = ∫ 0 t β S I t − s f s e τ + α s d s + ∫ t ∞ ϕ s − t f s e τ + α s ξ s − t d s − γ Q a t , d d t R t = α I t + γ Q c t + Q a t . (8) In Model (8), the spread of the disease relates to the number of contacts between susceptible and infected individuals. Therefore, we use Model (9), a simplified form of Model (8), to study the dynamics of the spreading process. d d t S t = − β S I t , d d t I t = β S I t − τ + α I t − ∫ 0 t β S I t − s f s e τ + α s d s − ∫ t ∞ f s ϕ s − t ξ s − t e τ + α t d s , d d t S S t = − 2 β m − 1 m S S t S I t S t , d d t S I t = β m − 1 m S S t S I t S t − β m − 1 m S I t S I t S t − β + τ + α S I t − ∫ 0 t β m − 1 m S S t − s S I t − s S t − s e ∫ t − s t β m − 1 m S I u S u + β + τ + α d u f s d s − ∫ t ∞ m S 0 ϕ s − t f s N ξ s − t e ∫ 0 t β m − 1 m S I u S u + β + τ + α d u d s . (9)

In this section, we determine the feasible solution region, calculate the basic reproduction number, and derive the final epidemic size of Model (9).

In this section, we apply the pairwise network model to determine the threshold of contact tracing rate and the effects of non-Markovian awareness process on epidemics theoretically.

To illustrate the accuracy of the analytical results, we compare the effects of different contact tracing rates and non-Markovian awareness processes by numerical simulations. The numerical simulations are performed with N = 1000, β = 0.1, α = 1/20, γ = 1/3, and m = 10. The initial conditions are [S](0) = 999, [I](0) = 1, [Q _c ](0) = [Q _a ](0) = [R](0) = 0, [ S S ] ( 0 ) = m N S 2 ( 0 ) and [SI](0) = m[S](0)[I](0)/N. According to [44], we assume the initial distribution of awareness process is ϕ ( s ) = [ Q a ] 0 δ ( s ) , where δ(s) is the Dirac delta function. Then, ∫ t ∞ ϕ ( s − t ) e − ( α + τ ) s f ( s ) ξ ( s − t ) d s = [ Q a ] 0 f ( t ) e − ( α + τ ) t . The awareness process of infection among infectious individuals is assumed to follow Gamma or Uniform distribution. Here, Gamma distribution is widely used in the epidemiology literature to approximate empirically observed transition periods [29]. For Gamma distribution Gamma(a, b), a is the shape parameter and b is the scale parameter. The probability density function is f ( x ) = x a − 1 e − x / b Γ ( a ) b a , when x > 0, where the function Γ(x) is Gamma function. The expected value is ab and the variance is ab ². Note that Gamma(1, λ) is Exponential distribution E(λ). Further, U(a, b) is a uniform distribution in the interval [a, b] with a ≥ 0 and a < b. The probability density function is f ( x ) = 1 b − a when x ∈ [a, b], and the expected value and the variance are (a+b)/2 and (b−a)²/12, respectively.

First, we focus on the effects of different contact tracing rates and awareness processes on the epidemic spreading.

Figure 2 shows the impacts of parameters on the basic reproduction number, R ₀. The contours of R ₀ as the function of the contact tracing rate and distribution parameters of awareness process are provided in Figure 2A, Figure 2B and Figure 2C. These simulations indicate that the basic reproduction number decreases as the contact tracing rate, τ, increases and the mean of awareness process decreases. Besides, the reduction in R ₀ that achieved by increasing τ is smaller when τ is larger. This means that it is relatively difficult to reduce the basic reproduction number by increasing the contact tracing rate when the contact tracing rate is larger. In Figure 2D, the basic reproduction numbers are compared when awareness processes have the same mean but different variances. We can observe that when the mean of awareness process is fixed, higher variance in the awareness process leads to smaller basic reproduction number and larger threshold of contact tracing rate. The simulation results in Figure 2 are consistent with the results of our theoretical analysis.

Figure 3 shows the impacts of parameters on the final epidemic size. In Figure 3A, Figure 3B and Figure 3C, the contours of the final epidemic size are depicted as the function of contact tracing rate and distribution parameters of awareness process. We find that the trends of these contours are similar, that is, the final epidemic size increases when the contact tracing rate increases and the mean awareness period decreases. Besides, the increase in the final epidemic size by increasing τ is greater when τ is larger. This indicates that it is relatively easy to increase the final epidemic size by increasing the contact tracing rate when the contact tracing rate is large. In Figure 3D, the awareness processes have the same mean but different variances. The results indicate that larger variance in awareness process contributes to larger final epidemic size when the tracing rate and the mean of awareness process are fixed.

The proportions of infection in the total population are shown in Figure 4 for different contact tracing rates and different distributions of the awareness process all having the same mean. Two observations can be obtained. First, larger contact tracing rate leads to a smaller and later peak of the proportion of infected individuals. Second, obvious differences exist in proportions of infected individuals despite the mean length of awareness periods is the same. In particular, when the contact tracing rate is fixed, the exponential distribution has the largest variance and leads to the lowest and latest peak in infection proportion.

Further, we are interested in the dynamics of contact behavior and awareness behavior of humans under different contact tracing rates and awareness processes. Here, we consider two non-Markovian distributions of awareness processes, U(0, 20) and Gamma(10, 1), which have the same mean but different variances. Note that the variance of U(0, 20) is larger than that of Gamma(10, 1).

The number of contacts between susceptible and infected individuals can be reduced by increasing contact tracing rates and the variance of awareness processes. In particular, Figure 5A and Figure 5B show that the number of S-I links has a lower and later peak when the contact tracing rate and the variance of awareness process are larger. The average number of contacts per susceptible individual with infectious individuals has similar characteristics, which is shown in Figure 5C and Figure 5D, respectively. The results indicate that a susceptible individual has less contacts with infectious individuals when the tracing rate is larger or the awareness process has a fixed mean but a larger variance.

Considering the variation of human awareness, we investigate the impacts of parameters on the proportions of aware individuals and traced individuals in Figure 6. Here, the proportion of aware individuals is the ratio of the cumulative number of aware individuals to the cumulative number of infected individuals, and the proportion of traced individuals is the ratio of the cumulative number of traced individuals to the cumulative number of infected individuals. We find that the proportion of aware individuals is larger when the awareness process has a fixed mean but a larger variance. However, the proportion of aware individuals is lower when the contact tracing is larger, which is possible because more infected individuals are traced in this case.

Finally, we investigate the performance of different contact tracing rates and awareness processes in containing epidemic spread in Table 1. Here, five characteristics of the epidemic, namely, the peak value of infected individuals, the total number of the infected individuals, the peak value of [SI], the total number of S-I links, and the peak value of [SI]/[S], are considered. The results indicate that larger reductions can be obtained by using larger contact tracing rates when the awareness process has a larger variance and a fixed mean, which is the combined effects of contact tracing and awareness process on the epidemic. Hence, enhancing contact tracing is more effective in reducing the number of infected individuals and the contacts between susceptible and infected individuals, when the awareness process has larger variance and a fixed mean.

The objective of this work is to determine the impacts of contact tracing strategy and awareness behavior among infectious individuals on the epidemic. Due to the complexity of human behavior, the non-Markovian awareness process is incorporated into the pairwise model in this work.

Through analyzing the SIR pairwise network model theoretically, we derive the basic reproduction number and an expression of final epidemic size. Moreover, we present the threshold for effective contact tracing rate, such that the basic reproduction number is less than one if the contact tracing rate is larger than the threshold. We prove that when the mean awareness period is fixed, a larger variance in awareness process contributes to a smaller basic reproduction number, a larger final epidemic size and a greater threshold of effective contact tracing rate. Further, two common non-Markovian distributions (Gamma and Uniform) are considered in numerical simulations to verify the analytic results.

We find that the enforcement of contact tracing has great impacts on the transmission of airborne diseases. Increasing the contact tracing rate can decrease the basic reproduction number, increase the final epidemic size, reduce the number of infections, thereby control the spread of disease. When the contact tracing rate is large, increasing contact tracing rate can greatly enlarge the final epidemic size. When the contact tracing rate is small, increasing tracing rate can reduce the basic reproduction number by a larger scale. To reduce the basic reproduction number to be less than one and finally curb the epidemic spread, an effective tracing rate that is larger than the threshold can be employed. Moreover, we find that susceptible individuals reduce their contacts with infectious individuals and more infectious individuals are traced when the contact tracing rate is larger. Hence, effective tracing strategies contribute to contain epidemic spread and reduce the contacts of susceptible individuals with infectious individuals, while the efficiencies of contact tracing vary.

The awareness on the state of infection among infected individuals is another important factor that affects the efficiency of disease control. When the mean awareness period is fixed, larger variance in the awareness process indicates that more infected individuals are aware of their infection in a shorter time. In this case, larger variance in the awareness process leads to smaller basic reproduction number, lower threshold of effective contact tracing rate, and larger treated (aware and traced) proportion of the infected population. Our study suggests that the public health authorities should improve the availability of rapid antigen testing and raise awareness of infectious status, so that more infected individuals are aware of their infected states and seek treatment in a short time. In particular, when the awareness process of infectious status has a larger variance and a fixed mean, contact tracing strategies have better performances in mitigating the spread of disease.

Our study incorporates contact tracing and the non-Markovian awareness process in a pairwise model. The analytical and simulation results help us understand the roles of contact tracing and awareness of humans in controlling the spread of disease. Since the economic costs of implementing contact tracing strategies are high especially when there are a large number of contacts with infected and treated people, public health authorities should also educate the population to maintain awareness about the infection. Rational allocation of resources for implementing contact tracing and awareness education can lead to rapid identification of infections and relieve the financial burden on the health authorities. In addition, our work can be extended by incorporating the non-Markovian infection process and recovery process into the epidemic model, which makes the model more realistic. Furthermore, the homogeneous network in our work is an approximation of the actual contact network. Hence, modelling heterogenous networks to investigate the role of contact tracing and awareness education in controlling the epidemic of infectious diseases will also be studied in our future work.