An SIS epidemic model with vaccination in a dynamical contact network of mobile individuals with heterogeneous spatial constraints

TL;DR

This paper models an SIS epidemic with vaccination in a dynamic network of moving individuals with heterogeneous interaction ranges, revealing how radius distribution affects disease spread and control strategies.

Contribution

It introduces a novel epidemic model incorporating heterogeneous interaction radii and mobility, deriving the basic reproduction number and analyzing vaccination effects in dynamic networks.

Findings

Heterogeneous radius distribution increases outbreak likelihood.

More heterogeneity leads to lower epidemic prevalence.

Optimal vaccination strategies are effective regardless of radius distribution.

Abstract

Network-based epidemic models have been extensively employed to understand the spread of infectious diseases, but have generally overlooked the fact that most realistic networks are dynamical rather than static. In this paper, we study a susceptible-infected-susceptible epidemic model with vaccination in a dynamical contact network of moving individuals, where we regard mobile individuals as random walkers that are allowed to perform long-range jumps. Different from previous studies of epidemics in a random walk network with a constant interaction radius, we consider the scenario where the individuals have a heterogeneous distribution of interaction radius $r$ and infected individuals are vaccinated with a probability depending on the interaction radius distribution. We derive the basic reproduction number $\mathcal{R}_0$ and explore the stability of disease-free and endemic equilibria of the model. Both theoretical and simulation results reveal that the distribution of individual interaction radius has significant effects on the basic reproduction number and the epidemic prevalence. In general, the disease will break out more readily in the population with a more heterogeneous radius distribution; however resulting in a lower epidemic prevalence. Interestingly, the results suggest that an optimal vaccination intervention for disease prevention and control is achievable regardless of the radius distribution. Furthermore, some interesting results on the structure of the underlying contact network are shown to have strong correlation with the epidemic dynamics. This study provides potential implications for developing efficient containment measures against infectious disease while considering the spatial constraints of moving individuals.

Equations

$$\begin{array}[]{ll}\displaystyle x_{i}(t+1)=x_{i}(t)+v\cos\eta_{i}(t),\\ \displaystyle y_{i}(t+1)=y_{i}(t)+v\sin\eta_{i}(t),\\ \displaystyle\eta_{i}(t+1)=\varepsilon_{i}(t+1),\\ \end{array}$$

$$d_{ij}(t)=d_{ji}(t)=\sqrt{\big{(}x_{i}(t)-x_{j}(t)\big{)}^{2}+\big{(}y_{i}(t)-y_{j}(t)\big{)}^{2}},\quad i,j=1,2,\dots,N.$$

$$\theta (r)={\theta }_0\frac{P(r)r^{\alpha }}{{\sum }_{j=1}^{m}P(r_j){r_j}^{\alpha }},$$

$$\begin{array}[]{ll}\displaystyle\frac{ds(r,t)}{dt}=\lambda i(r,t)+\varphi v(r,t)-s(r,t)\theta(r)-s(r,t)[1-(1-\beta)^{k_{r}^{inf}}],\\ \displaystyle\frac{di(r,t)}{dt}=-\lambda i(r,t)+s(r,t)[1-(1-\beta)^{k_{r}^{inf}}],\\ \displaystyle\frac{dv(r,t)}{dt}=s(r,t)\theta(r)-\varphi v(r,t),\quad r=r_{1},r_{2},\dots,r_{m}.\end{array}$$

$$\begin{array}[]{ll}\displaystyle\frac{ds(r,t)}{dt}=\lambda i(r,t)-s(r,t)\beta\rho\pi r^{2}\sum_{j=1}^{m}P(r_{j})i(r_{j},t)+\varphi v(r,t)-s(r,t)\theta(r)\\ \displaystyle\frac{di(r,t)}{dt}=-\lambda i(r,t)+s(r,t)\beta\rho\pi r^{2}\sum_{j=1}^{m}P(r_{j})i(r_{j},t)\\ \displaystyle\frac{dv(r,t)}{dt}=s(r,t)\theta(r)-\varphi v(r,t),\quad r=r_{1},r_{2},\dots,r_{m}.\end{array}$$

$$\begin{array}[]{ll}s(r,0)+i(r,0)+v(r,0)=1,\\ 0\leq s(r,0)\leq 1,\\ 0\leq i(r,0)\leq 1,\\ 0\leq v(r,0)\leq 1,\quad r=r_{1},r_{2},\dots,r_{m}.\end{array}$$

$$0\le s(r,t),i(r,t),v(r,t)\le 1,s(r,t)+i(r,t)+v(r,t)=1$$

$$t_{0}=inf\{t\geq 0\big{|}i(r_{0},t)=0\},$$

$$\frac{di(r_0,t_0)}{dt}=s(r_0,t_0)\beta \rho \pi r_0^2\sum _{j=1}^{m}P(r_j)i(r_j,t_0)<0.$$

$$\frac{ds(r_0,t_1)}{dt}=\lambda i(r_0,t_1)+\phi v(r_0,t_1)<0,$$

$$\frac{dv(r_0,t_2)}{dt}=s(r_0,t_2)\theta (r_0)<0.$$

$$t_{0}=inf\{t\geq 0\big{|}s(r_{0},t)=0\},$$

$$\frac{ds(r_0,t_0)}{dt}=\lambda i(r_0,t_0)+\phi v(r_0,t_0)<0,$$

$$\begin{array}[]{ll}\displaystyle\frac{ds(r,t)}{dt}=&\lambda i(r,t)-s(r,t)\beta\rho\pi r^{2}\sum_{j=1}^{m}P(r_{j})i(r_{j},t)-s(r,t)\theta(r)\\ &+\varphi\big{[}1-s(r,t)-i(r,t)\big{]},\\ \displaystyle\frac{di(r,t)}{dt}=&-\lambda i(r,t)+s(r,t)\beta\rho\pi r^{2}\sum_{j=1}^{m}P(r_{j})i(r_{j},t),\quad r=r_{1},r_{2},\dots,r_{m}.\end{array}$$

$$\begin{array}[]{ll}\Gamma=&\Bigg{\{}\Big{(}s(r_{1},t),i(r_{1},t),\dots,s(r_{m},t),i(r_{m},t)\Big{)}\in\mathbb{R}^{2m}\Big{|}0\leq s(r,t)\leq 1,\\ &0\leq i(r,t)\leq 1,0\leq s(r,t)+i(r,t)\leq 1,r=r_{1},r_{2},\dots,r_{m}\Bigg{\}},\end{array}$$

$$E_{0}=\Big{(}\frac{\varphi}{\varphi+\theta(r_{1})},0,\frac{\varphi}{\varphi+\theta(r_{2})},0,\dots,\frac{\varphi}{\varphi+\theta(r_{m})},0\Big{)},$$

$$\displaystyle i_{r}^{*}=\frac{{\beta\varphi\rho\pi r^{2}}{\sum_{j=1}^{m}P(r_{j})i_{r_{j}}^{*}}}{\lambda\big{(}\varphi+\theta(r)\big{)}+\beta\varphi\rho\pi r^{2}{\sum_{j=1}^{m}P(r_{j})i_{r_{j}}^{*}}},$$

$$I(t)=\sum _{j=1}^{m}P(r_j)i(r_j,t),I^{\ast }=\sum _{j=1}^{m}P(r_j)i_{r_j}^{\ast }.$$

$$\displaystyle I^{*}=\sum_{j=1}^{m}P(r_{j}){\frac{\beta\varphi\rho\pi r_{j}^{2}I^{*}}{\lambda\big{(}\varphi+\theta(r_{j})\big{)}+\beta\varphi\rho\pi r_{j}^{2}I^{*}}},$$

$$\displaystyle F(I^{*})=\sum_{j=1}^{m}P(r_{j}){\frac{\beta\varphi\rho\pi r_{j}^{2}I^{*}}{\lambda\big{(}\varphi+\theta(r_{j})\big{)}+\beta\varphi\rho\pi r_{j}^{2}I^{*}}}-I^{*},$$

$$\displaystyle\frac{dF(I^{*})}{dI^{*}}=\sum_{j=1}^{m}P(r_{j}){\frac{\beta\varphi\rho\pi r_{j}^{2}\big{(}\lambda\varphi+\lambda\theta(r_{j})\big{)}}{\big{(}\beta\varphi\rho\pi r_{j}^{2}I^{*}+\lambda\varphi+\lambda\theta(r_{j})\big{)}^{2}}}-1,$$

$$\displaystyle\frac{d^{2}F(I^{*})}{d{I^{*}}^{2}}=\sum_{j=1}^{m}P(r_{j}){\frac{-2(\beta\varphi\rho\pi r_{j}^{2})^{2}\big{(}\lambda\varphi+\lambda\theta(r_{j})\big{)}}{\big{(}\beta\varphi\rho\pi r_{j}^{2}I^{*}+\lambda\varphi+\lambda\theta(r_{j})\big{)}^{3}}}<0.$$

$$\displaystyle F(0)=0,F(1)=\sum_{j=1}^{m}P(r_{j}){\frac{\beta\varphi\rho\pi r_{j}^{2}}{\lambda\big{(}\varphi+\theta(r_{j})\big{)}+\beta\varphi\rho\pi r_{j}^{2}}}-1<0,$$

$$\displaystyle\frac{dF(I^{*})}{dI^{*}}\Bigg{|}_{I^{*}=0}=\sum_{j=1}^{m}P(r_{j}){\frac{\beta\varphi\rho\pi r_{j}^{2}}{\lambda(\varphi+\theta(r_{j}))}}-1>0.$$

$${\mathcal{R}}_0=\frac{\beta \phi \rho \pi }{\lambda }\sum _{j=1}^{m}P(r_j)\frac{r_j^2}{\phi +\theta (r_j)}.$$

$$E^{*}=\Big{(}s^{*}_{r_{1}},i^{*}_{r_{1}},s^{*}_{r_{2}},i^{*}_{r_{2}},\dots,s^{*}_{r_{m}},i^{*}_{r_{m}}\Big{)},$$

$$R_0=\frac{\beta \phi \rho \pi }{\lambda (\phi +{\theta }_0)}\sum _{j=1}^{m}P(r_j)r_j^2,$$

$$J_{E_0}={\left(\begin{array}{ccccc}{A}_1& B_{12}& B_{13}& \dots & B_{1m}\\ B_{21}& A_2& B_{23}& \dots & B_{2m}\\ ⋮& ⋮& ⋮& & ⋮\\ B_{m1}& B_{m2}& B_{m3}& \dots & A_m\end{array}\right)}_{2m\times 2m},$$

$$A_{j}=\begin{pmatrix}-\Big{(}\theta(r_{j})+\varphi\Big{)}&\lambda-\varphi-\frac{\varphi}{\varphi+\theta(r_{j})}\beta\rho\pi r_{j}^{2}P(r_{j})\\ 0&-\lambda+\frac{\varphi}{\varphi+\theta(r_{j})}\beta\rho\pi r_{j}^{2}P(r_{j})\end{pmatrix},\quad j=1,2,\dots,m,$$

$$B_{ij}=\left(\begin{array}{cc}0& -\frac{\phi }{\phi +\theta (r_i)}\beta \rho \pi r_i^{2}P(r_j)\\ 0& \frac{\phi }{\phi +\theta (r_i)}\beta \rho \pi r_i^{2}P(r_j)\end{array}\right),i,j=1,2,\dots ,m.$$

Full text

††journal: Elsevier

1 Introduction

As yet, a powerful technique to study infectious diseases spreading in populations has been to build a mathematical compartment model [1, 2, 3, 4] through which we can grasp how the disease evolves in time and provide implications for devising effective disease control measures. Two classical examples, among others, are the susceptible-infected-susceptible (SIS) model in which individuals can be reinfected after recovery and the susceptible-infected-recovered (SIR) model in which individuals gain lifelong immunity after infection [1, 2, 3, 4]. Though very useful for dynamical analysis, these compartment models nevertheless have heavily hinged upon the homogeneous mixing (HM) assumption that all the individuals have the same probability to contact each other [1, 2, 3, 4]. However, the contact pattern among individuals in real populations is heterogeneous such that the number of contacts or interactions varies widely among individuals [5, 6]. In recent decades, complex networks [5, 7, 8, 9] have become a repeatedly used paradigm to study the spread of infectious diseases since the network-based epidemiological models [10, 11, 12, 13, 14, 15, 16, 17, 18] have transcended the limitation of HM assumption in compartment models. The underlying network structures have been found to entail significant impacts on epidemic spreading dynamics [11] and other dynamical processes taking place over networks [19]. As a well-known example, Pastor-Satorras and Vespignani [20] discovered that the epidemic threshold for the SIS epidemic model vanishes in the thermodynamic limit in a static scale-free (SF) network [5] which obeys a power-law degree distribution $P(k)\propto k^{-\gamma}$ with an exponent $2<\gamma\leq 3$. It is worth remarking that a more mathematically rigorous analysis showed that the critical value is also zero for contact processes on random networks with power-law degree distributions for any value of power $\gamma>3$ [21].

A large body of investigations on network epidemiology have focused on static networks, where the topological structure of network is fixed in time during the process of disease transmission [22, 23, 24, 25, 26, 27, 28, 29, 30, 31], thus ignoring the impact of motions of individuals. In reality, individuals in a population often move around during the spread of infectious diseases, leading to structural changes of the underlying contact network that mimics the population [32, 33]. Recently, the epidemic spreading on random walk networks [34, 35, 36, 37, 38, 39, 40, 41] has been extensively investigated to understand the effects of individuals’ motion on the epidemic dynamics. For instance, Frasca et al. [35] proposed a dynamical network model with mobile individuals who are allowed to perform both local and long-distance motions. In their model, mobile individuals are modeled as random walkers who are only able to interact with others falling within a given interaction radius apart from them. Based on a similar dynamical network model, Buscarino et al. [36] argued that the homogeneous mixing approximation is appropriate only when the velocity of individuals’ movement is large enough.

On the one hand, most infectious disease models presented in the literature have largely neglected the influence of spatial distances between individuals. In fact, many realistic networks such as the mobile phone communication network, social contact network and the power grid are often embedded in a Euclidean geographical space [42] and the interactions among individuals usually depend on their spatial distances and geographical information [43]. Generally speaking, the living space and the sphere of activity of individuals are constrained in terms of spatial distances [44, 45]. Moreover, because of individual diversity, the activity ability, activity range, contact number and the geographic location [46] differ from individual to individual. Typically the limitation of the individuals on the spatial distance will lead to the localization of the contact pattern between individuals, which will affect the transmission of infectious diseases throughout the contact network [44, 45]. Therefore, it is natural to study the infectious diseases model in contact networks with geographical properties, such as the embedded lattice [43] and the spatially embedded networks [47]. Xu $et$ $al.$ presented an SIS epidemic model in a lattice-embedded scale-free network and investigated how the geographical structure affects the dynamical process of epidemic spreading [48]. As a further step, Xu $et$ $al.$ also considered the standard SIS model on a random growing network to study the integrated effects of preference and geography on epidemic spreading [49]. In such spatially embedded networks, the individuals’ interaction radius is generally assumed to be primarily determined by their respective degrees, that is, the larger degree, the larger interaction radius [43, 47, 48, 49]. It is worth remarking that most previous mentioned works based on random walk networks [34, 35, 36, 37] have simply assumed that all individuals have the same interaction radius in order to better include other factors such as the velocity and the direction of motion, as well as the population density. However, the interaction radius of individuals in realistic populations or networks are usually heterogeneous [50]. For example, individuals with poor personal hygiene are prone to have a larger radius of contacting infectious sources. In a wireless sensor communication network, sensors with different power have different communication radii [51]. Most of the aforementioned works concentrated on the threshold analysis of the model under consideration. However, there is little (if any) work devoted to stability analysis on epidemic models in complex networks with spatial or geographical constraints in the literature. In this paper, we consider an epidemic model with heterogeneous interaction radius of individuals, based on which we derive the basic reproduction number, analyze the equilibria stability, prove the model persistence, and investigate the effects of spatial constraints in individuals’ mobilities and vaccination intervention on the epidemic spreading and on the network structure.

On the other hand, vaccination is one of the most effective policies for preventing the transmission of infectious diseases [52, 53] and up to now there have been a large number of studies on various vaccination strategies for epidemic models in complex networks [54, 55, 56, 57, 58, 59]. Traditional vaccination methods include random and targeted immunization strategies [60]. It is argued that random immunization strategy is insufficient for networks with broad degree distributions, whereas targeted vaccination is to immunize high-degree nodes and has a much higher effectiveness than random vaccination in SF networks [60] and small-world networks [61]. However, it is difficult to implement the targeted vaccination strategy since it requires full knowledge of the degree of each node in the network. To overcome this shortcoming, several vaccination strategies based on local information have been proposed, such as acquaintance immunization [62]. All these mentioned vaccination strategies are based on the degree distribution of individuals of the network. Rather than relying on node degree, in this paper we adopt individuals’ interaction radius to characterize the individual heterogeneity and propose a new vaccination strategy that depends on the interaction radius of each individual.

The contributions of this paper are as follows. An SIS epidemic model with interaction radius-dependent vaccination is proposed to probe the potential effects of heterogeneous spatial constraints of individuals on epidemic spreading in a dynamic contact network of moving individuals. The epidemic dynamics is described by a set of ordinary differential equations. The explicit mathematical expression of the basic reproduction number is derived and the dynamical properties of both the disease-free equilibrium and the endemic equilibrium of the model system are presented. Numerical calculation and stochastic simulation show good agreement, indicating that our model can well describe the dynamical process of disease transmission on dynamical networks of mobile individuals. The effects of different radius distributions on epidemic dynamics and network structure have been examined. Our results imply that the optimal vaccination intervention is realizable for disease prevention and control.

This paper is outlined as follows. In section 2, we describe the construction of a random walk network in which individuals are allowed to perform long-distance jumps with a probability and then present the SIS epidemic model with vaccination that depends on individuals’ interaction radius. In section 3 we derive the basic reproduction number $\mathcal{R}_{0}$. We give the stability analysis of equilibria in Section 4 and discuss the persistence of the model in Section 5. In section 6 we give the simulation results and discussion. In section 7, we conclude the paper.

2 Model description

2.1 Dynamic contact network of mobile individuals

There are a number of works that rely on random walk network to inspect the effects of moving agents on epidemic spreading [34, 35, 36, 39, 50]. In a similar framework, we consider $N$ individuals who are initially randomly distributed in a two-dimensional space $\Omega=\{{(x,y)}\in\mathbb{R}^{2}:0\leq x\leq D,0\leq y\leq D\}$, with periodic boundary conditions, as illustrated in Fig. 1.

For convenience, we denote $\Lambda_{i}(t)=(x_{i}(t),y_{i}(t))$ as the position of the individual $i$ $(i=1,2,\dots,N)$ in the planar space with its moving velocity $v_{i}(t)=(v\cos\eta_{i}(t),v\sin\eta_{i}(t))$ and moving direction $\eta_{i}(t)$ at time $t$, where $v$ is the modulus of the agent velocity, which is the same for all individuals. Then, the motion of individual $i$ can be described as follows:

[TABLE]

where $\varepsilon_{i}(t)$ is a random variable obeying the uniform distribution between the interval $[-\pi,\pi]$. In addition, to include the probability that individuals can move with time scales much shorter than those related to disease, we consider the case where infected individuals may perform long-distance jumps. We define a parameter $p_{\rm{jump}}$ that quantifies the probability for an individual to perform a long-distance jump. Each individual can jump to any position (i.e. long-distance jump) inside the planar space with the probability $p_{\rm{jump}}$, that is, $p_{\rm{jump}}$ denotes the probability of an individual jumping to a random position in the space $\Omega$, similar to the case of other works [34, 35, 36, 50]. Each individual evolves following $v_{i}(t)=(v\cos\eta_{i}(t),v\sin\eta_{i}(t))$ with probability $1-p_{\rm{jump}}$ or performs a random jump with probability $p_{\rm{jump}}$. In what follows, the model is investigated as a function of the parameter $p_{\rm{jump}}$. At time $t$, the Euclidean distance between individual $i$ and $j$ is defined as

[TABLE]

In the present work, we consider that each individual’s behavior is constrained by their respective spatial distance, which is characterized by the interaction radius $r$ in the following. In order to take account of individual heterogeneity, we consider an individual $i$ has its own interaction radius $r^{i}$, i=1,2,…,N. Here, the interaction radius is the effective interaction distance on the Euclidean plane, denoting the characteristic radius of the circular region within which individuals can get infections from others. In our model, we assume there are $m$ different values of interaction radius which obey a preassigned probability distribution $P(r_{j})$, $j=1,\dots,m$, where $P(r_{j})$ denotes the proportion of nodes with interaction radius $r_{j}$. Namely, the interaction radius of each individual is given once and for all. At any time $t$, individual $i$ can only be able to interact with other individuals $j(\neq i)$ that fall within the circle defined by individual $i$’s position $\Lambda_{i}(t)$ and its interaction radius $r^{i}$. In the context of disease transmission, individual $i$ can be infected by any infected individuals who are located within the circle defined by the location $\Lambda_{i}$ and interaction radius $r^{i}$ of individual $i$. In this sense, the interaction radius can be seen as the “susceptibility” radius [50]. As demonstrated in Fig. 1, it is possible for the disease to spread from the individual $1$ to individual $6$, whereas it is impossible to spread from individual $9$ to individual $5$ since the individual $9$ does not enter the realm of individual $5$. The definition of interaction radius forms a dynamical directed contact network [63, 64], as illustrated in Fig. 2. All the individuals comprise the nodes of the contact network, in which the contacts are defined asymmetrically such that a node $j$ is regarded as a effective contact or a neighbor (which is capable of disease transmission) of node $i$ only if node $j$ is located in the circular realm of node $i$, but the converse is not necessarily true. In the layout given in Fig. 1, for example, both the nodes $3$ and $7$ are neighbors of node $2$, while the node $2$ is neither the neighbor of node $3$ nor the neighbor of node $7$. Nodes $4$ and $8$ are mutual effective contacts to each other since either of them is positioned within the interaction radius of the other. Note that when all the individuals share an identical interaction radius, i.e., the values of interaction radius follow a delta distribution. The in-degree and out-degree of nodes of the directed network [63, 64] can be defined as follows. The in-degree of node $i$ at time $t$ depends on the number of other nodes that fall into the interaction radius $r^{i}$ of node $i$. That is, the in-degree of node $i$ is defined as the number of effective contacts in the realm of node $i$. More clearly, the in-degree of node $i$ is the number of nodes $j$ who satisfy $d_{ij}(t)\leq r^{i},\forall j\neq i$. Similarly, the out-degree of node $i$ is defined as the number of nodes whose effective contacts include node $i$, that is, the number of nodes $j$ who meet $d_{ij}(t)\leq r^{j},\forall j\neq i$.

2.2 SIS epidemic model with vaccination

In this paper, we consider an SIS epidemic model with vaccination, where the vaccination of susceptible individuals depends on their interaction radius. The choice of such a vaccination strategy is motivated by the following consideration. We observe that a lot of studies on vaccination strategies take the node degree as the characteristic index of individual heterogeneity [22, 23, 60, 62, 65, 66, 67, 68, 69]. As an alternative measurement to quantify the diversity among individuals, interaction radius can characterize the range and ability of individuals’ activity [50]. In this regard, we consider different vaccination rate for susceptible individuals according to their interaction radius.

In the model, there are a total number $N$ of individuals, each of which may have only one of the three possible states: susceptible (S), infected (I) and vaccinated (V). We denote $N_{S}(t)$, $N_{I}(t)$ and $N_{V}(t)$ as, respectively, the number of susceptible, infected and vaccinated individuals at time $t$. Obviously, we have the total number of individuals $N_{S}(t)+N_{I}(t)+N_{V}(t)=N$, which remains constant over time.

The epidemic spreading process follows the transmission rule as shown in Fig. 3, where the individuals of states S, I and V are all subscripted with an interaction radius $r$. At each time step, each susceptible individual can be infected by an infected neighbor with the transmission rate $\beta$ and be vaccinated with the vaccination rate $\theta(r)$ that depends on the interaction radius of the susceptible individual. Each infected individual recovers at the recovery rate $\lambda$, and each vaccinated individual returns to being susceptible again with the resusceptibility rate $\varphi$ after the vaccine wears off.

Here, we consider the vaccination rate $\theta(r)$ in the form

[TABLE]

where $\alpha\in\mathbb{R}$ is a tunable parameter and $0<\theta_{0}<1$ is a constant. Equation (1) implies that the vaccination of an individual depends on the individual’s effective interaction radius $r$ and the probability $P(r)$ that the individual’s interaction radius is $r$. When $\alpha=0$, it means that the vaccination rate depends only on the probability distribution $P(r)$ of the interaction radius. When $\alpha\neq 0$, the vaccination rate is $\theta(r)\propto P(r)r^{\alpha}$, similar to the preferential attachment hypothesis in the growing network model [70]. The case of $\alpha>0$ means that the larger the value of the radius $r$ and its probability $P(r)$, the higher the vaccination rate; while the case of $\alpha<0$ indicates that the larger the value of $r$ and the smaller the probability $P(r)$, the smaller the vaccination rate. In the context of disease transmission, the case of $\alpha>0$ accounts for the scenario in which individuals who are active in social contacts have a higher risk of infections and should be vaccinated with high priority. On the contrary, the case of $\alpha<0$ considers the situation that individuals with close contacts should be protected preferentially.

Aiming to understand the system behavior of the epidemic spreading over the dynamical contact network of mobile individuals, we derive a mean-field model based on the HM assumption that the population mixes at random, i.e., all the individuals have the same probability to contact the other ones [1, 2, 3, 4]. It has been argued that if the individuals’ mobility is extremely high, namely, if $p_{\rm{jump}}\rightarrow 1$ or $v\simeq D$ is satisfied, the underlying contact network is degenerated into an averaged one that can be well approximated by the HM hypothesis [34, 35, 36, 50]. In our model, the parameter $p_{\rm{jump}}$ quantifies the probability for an individual to perform a long-distance jump to a random position in the planar space. When $p_{\rm{jump}}=1$, all the individuals jump to a random position in the two-dimensional space independently. In this case, the spatial correlations in the disease states are destroyed by the agent motion [71]. Therefore, the case of $p_{\rm{jump}}=1$ is equivalent to the case where all the individuals have the same chance to encounter the others. In a word, the case of $p_{\rm{jump}}=1$ in our model can be well approximated by the HM assumption.

Under the HM hypothesis, we have the density of individuals $\rho=N/{D^{2}}$ which can be fixed by fixing the value of $D$. Taking into account the heterogeneous interaction radius of individuals, we denote the number of susceptible, infected and vaccinated individuals with effective interaction (susceptibility) radius $r$ by $N_{S}(r,t)$, $N_{I}(r,t)$ and $N_{V}(r,t)$, respectively. It is straightforward to get the number of individuals with effective interaction radius $r$ as $N_{r}=N_{S}(r,t)+N_{I}(r,t)+N_{V}(r,t)$ which is a constant only depending on the probability distribution $P(r)$. We further define $s(r,t)=N_{S}(r,t)/N_{r}$, $i(r,t)=N_{I}(r,t)/N_{r}$ and $v(r,t)=N_{V}(r,t)/N_{r}$ as the fraction (or relative density) of individuals with effective interaction radius $r$ at time $t$, respectively in the susceptible, infected and vaccinated states. Obviously, $s(r,t)+i(r,t)+v(r,t)=1$. The equations for the dynamical system of mobile individuals in different disease states read

[TABLE]

The first equation refers to creation of susceptibles because of recovery of infecteds (first term), increase of susceptibles due to resusceptibility of vaccinateds (second term), loss of susceptibles due to vaccination (third term) and decrease of susceptibles caused by infection that is proportional to $s(r,t)$ times a contagion probability $p_{\rm{cont}}$. The contagion probability is given by $p_{\rm{cont}}=1-\bar{p}_{\rm{cont}}$, where $\bar{p}_{\rm{cont}}$ is the probability of not being infected. Since $\bar{p}_{\rm{cont}}$ is the probability that an individual with radius $r$ is not infected by any of its infected neighbors at time $t$, we have $\bar{p}_{\rm{cont}}=(1-\beta)^{k_{r}^{inf}}$, where $k_{r}^{inf}$ is the number of its infected neighbors. Based on the mean-field approximation and the HM hypothesis, we have $k_{r}^{inf}=\rho\pi r^{2}\sum_{j=1}^{m}P(r_{j})i(r_{j},t)$ [50]. Therefore, the fraction of susceptible individuals with radius $r$ that enter the infected compartment at time $t+1$ is $s(r,t)[1-(1-\beta)^{k_{r}^{inf}}]$. The second equation indicates that the decrease of infecteds is proportional to the fraction of infecteds which get recovered, i.e., to $\lambda i(r,t)$, and that the increase of infecteds is proportional to the density of susceptibles contracting the disease. The third equation is derived by considering the vaccination of susceptibles with rate $\theta(r)$ and the relapse into susceptibility with rate $\varphi$ for vaccinated individuals as the vaccine wears off.

As the transmission rate $\beta$ is small enough, we can make the approximation $[1-(1-\beta)^{k_{r}^{inf}}]\approx\beta{k_{r}^{inf}}$. In this paper, the epidemiological parameters are set to be small since the choice of parameter values only affects the time scale of the disease propagation without influencing the generality of the results [37, 72, 73]. Following such an approximation, the model system of Eqs. (2) can be simplified to

[TABLE]

The initial conditions of the model system (3) take the form

[TABLE]

Combined with the initial conditions (4) and the preassigned interaction radius distribution $P(r)$, the system (3) determines the epidemic dynamics on the spatial contact network of mobile individuals. In order to present the dynamical properties of the system, in what follows we provide a lemma on the positiveness and boundedness of solutions to system (3).

Lemma 1**.**

Let \big{(}s(r,t),i(r,t),v(r_{,}t)\big{)},r=r_{1},r_{2},\dots,r_{m} be the solutions to system (3) with the initial conditions given by (4), then it follows that

[TABLE]

for any $r=r_{1},r_{2},\dots,r_{m}$ and $t\geq 0$.

Proof 1**.**

Firstly, we verify $i(r,t)\geq 0$ for any $r=r_{1},r_{2},\dots,r_{m}$. By the way of contradiction, because $i(r,0)\geq 0$, we assume that there exist some $r_{0}\in\{r_{1},r_{2},\dots,r_{m}\}$ and $t\geq 0$ such that $i(r_{0},t)=0$. Let

[TABLE]

then $i(r_{0},t_{0})=0$, $\frac{di(r_{0},t_{0})}{dt}<0$ and $i(r_{0},t)>0$ for any $t\in[0,t_{0})$. It follows from the second equation of system (3) that

[TABLE]

This indicates that $s(r_{0},t_{0})<0$. Since $s(r_{0},0)\geq 0$, there exists a $t_{1}<t_{0}$ such that $s(r_{0},t_{1})=0$, $\frac{ds(r_{0},t_{1})}{dt}<0$ and $s(r_{0},t)>0$ for any $t\in[0,t_{1})$.

By the first equation of system (3), we have

[TABLE]

which implies that $v(r_{0},t_{1})<0$, since $i(r_{0},t_{1})>0$.

Similarly, there exists a $t_{2}<t_{1}$ such that $v(r_{0},t_{2})=0$, $\frac{dv(r_{0},t_{2})}{dt}<0$ and $v(r_{0},t)>0$ for any $t\in[0,t_{2})$. Substituting $v(r_{0},t_{2})=0$ into the last equation of system (3) yields

[TABLE]

This means $s(r_{0},t_{2})<0$ which leads to a contradiction with $s(r_{0},t)>0$ for any $t\in[0,t_{1})$. Hence, $i(r,t)\geq 0$ for any $r=r_{1},r_{2},\dots,r_{m}$ and any $t\geq 0$.

In a similar way, using the way of contradiction starting from the third equation of system (3), we can easily show that $v(r,t)\geq 0$ for any $r=r_{1},r_{2},\dots,r_{m}$ and any $t\geq 0$.

Now we assume there exist some $r_{0}\in\{r_{1},r_{2},\dots,r_{m}\}$ and $t\geq 0$ such that $s(r_{0},t)=0$. Let

[TABLE]

then $s(r_{0},t_{0})=0$, $\frac{ds(r_{0},t_{0})}{dt}<0$ and $s(r_{0},t)>0$ for any $t\in[0,t_{0})$. It follows from the first equation of system (3) that

[TABLE]

which is a contraction since $i(r_{0},t_{0})\geq 0$ and $v(r_{0},t_{0})\geq 0$. This contradiction indicates that $s(r,t)\geq 0$ for any $r=r_{1},r_{2},\dots,r_{m}$ and any $t\geq 0$.

Therefore, it is straightforward to obtain $s(r,t)\leq 1$, $i(r,t)\leq 1$, and $v(r,t)\leq 1$ because $s(r,t)+i(r,t)+v(r,t)=1$. The proof of Lemma 1 is completed. ∎

3 Equilibria and basic reproduction number

In this section, we will derive the basic reproduction number $\mathcal{R}_{0}$ by examining the existence and uniqueness of the endemic equilibrium (EE) of our model. In mathematical epidemiology, the basic reproduction number is an important threshold indicator that determines whether the disease breaks out or dies out. It is defined as the average number of new infections caused by an infected individual during its infectious period when appearing in a completely susceptible population [1, 2, 3, 4]. Generally, if $\mathcal{R}_{0}>1$ then the disease will break out resulting in an endemic state; otherwise if $\mathcal{R}_{0}<1$ the disease will become extinct eventually [1, 2, 3, 4].

Based on the normalization condition $s(r,t)+i(r,t)+v(r,t)=1$, the model (3) can be reduced to

[TABLE]

According to Lemma 1, the feasible region for system (5) is given by

[TABLE]

which is positively invariant with regard to system (5).

Obviously, system (5) admits a unique disease-free equilibrium (DFE)

[TABLE]

on the boundary $\partial\Gamma$ of the invariant set $\Gamma$.

By letting the right-hand side of (5) be zero, we have a stationary solution of system (5) in the limit of $t\rightarrow\infty$ as

[TABLE]

where $i_{r}^{*}=\lim\limits_{t\rightarrow\infty}i(r,t)$, $r=r_{1},r_{2},\dots,r_{m}$.

Denote by $I(t)\in[0,1]$ the fraction of infected individuals among the total $N$ individuals and by $I^{*}$ the stationary value of $I(t)$ as $t\rightarrow\infty$. By this definition, we have

[TABLE]

Note that $0\leq i(r,t)\leq 1$ for any $r=r_{1},r_{2},\dots,r_{m}$ and $t\geq 0$ according to Lemma 1, thus we have $0\leq i_{r}^{*}\leq 1$ for any $r=r_{1},r_{2},\dots,r_{m}$, and hence $0\leq I^{*}\leq 1$.

Combining Eqs. (8) and (9) gives rise to a self-consistency equation

[TABLE]

which implies a trivial solution $I^{*}=0$. Now we give the conditions about the existence and uniqueness of the nontrivial positive solution $I^{*}>0$. To this aim, define

[TABLE]

then we have

[TABLE]

and

[TABLE]

That is, the continuous function $F(I^{*})$ is convex upward in the interval $[0,1]$. In addition, since

[TABLE]

then it follows from the continuity of the function that the necessary and sufficient condition for the existence and uniqueness of the positive solution $0<I^{*}<1$ to Eq. (10) should be

[TABLE]

This inequality determines the basic reproduction number $\mathcal{R}_{0}$ of our model as follows:

[TABLE]

In other words, when $\mathcal{R}_{0}>1$, there exists a unique positive solution $0<I^{*}<1$ to Eq. (10) in addition to the trivial solution $I^{*}=0$; otherwise, the trivial solution $I^{*}=0$ is the only solution to Eq. (10). Furthermore, given that $0<I^{*}<1$, it follows from Eq. (8) that $0<i^{*}_{r}<1$ for any $r=r_{1},r_{2},\dots,r_{m}$. Therefore, if $\mathcal{R}_{0}>1$, then the model system (5) has a unique EE point $E^{*}$ given as

[TABLE]

where $0<s^{*}_{r}<1$, $0<i^{*}_{r}<1$ for any $r=r_{1},r_{2},\dots,r_{m}$. We summarize the above statements in the following theorem.

Theorem 2**.**

Consider the system (5) and the basic reproduction number $\mathcal{R}_{0}$ defined by Eq. (11), then there always exists a DFE point $E_{0}$ as given in Eq. (7). Moreover, if and only if $\mathcal{R}_{0}>1$, the system (5) as a unique EE point $E^{*}$ as given in Eq. (12). Therefore, the epidemic breaks out when $\mathcal{R}_{0}>1$; otherwise, the disease dies out eventually.

Remark 1**.**

In fact, the basic reproduction number $\mathcal{R}_{0}$ can be interpreted in the epidemiological perspective as follows. Consider a healthy population without any infected seed, then according to the model definition and Eq. (7), there are only susceptible individuals and vaccinated ones in the steady state, where the stationary number of susceptible individuals in each compartment with effective interaction radius $r=r_{1},r_{2},\dots,r_{m}$ is given by $N_{S}^{*}(r)=N\frac{\varphi}{\varphi+\theta(r)}$ with the probability distribution $P(r)$. When an infected seed is introduced in the population, the infected individual contacts a susceptible individual with interaction radius $r$ with probability $\pi r^{2}/D^{2}$ and transmits the disease to the susceptible at rate $\beta$. Therefore, on average, the infected seed will create $\frac{\beta}{\lambda}\sum_{j=1}^{m}P(r_{j})N\frac{\varphi}{\varphi+\theta(r_{j})}\frac{\pi r_{j}^{2}}{D^{2}}$ new infections during its entire infectious period $\tau=1/\lambda$. Using the definition of $\rho$ gives $\mathcal{R}_{0}$.

Remark 2**.**

In the extreme case of $\theta(r)=\theta_{0}$ where a random vaccination scheme is adopted, the basic reproduction number $\mathcal{R}_{0}$ is given by

[TABLE]

indicating that $\mathcal{R}_{0}$ is proportional to the transmission rate $\beta$, the population density $\rho$ and the second moment $\langle r^{2}\rangle=\sum_{j=1}^{m}P(r_{j})r_{j}^{2}$ of the radius distribution. In particular, if $\theta_{0}=0$, then there is no vaccination and our model reduces to SIS model, where we get $R_{0}={\beta\rho\pi}\langle r^{2}\rangle/{\lambda}$, reproducing the result obtained in [50].

4 Stability analysis

In this section, we will study the local and global dynamics of DFE point $E_{0}$ given by (7) and EE point $E^{*}$ given by (12) of the model system (5). We present all the results on the dynamical behavior of the equilibria in the following theorems.

4.1 Stability of DFE

Theorem 3**.**

Consider the model system (5), the following two conclusions hold.

• (1)

If $\mathcal{R}_{0}<1$, then the DFE point $E_{0}$ is locally asymptotically stable.

• (2)

If $\mathcal{R}_{0}>1$, then the DFE point $E_{0}$ is unstable.

Proof 2**.**

To determine the local stability of DFE point $E_{0}$, we consider the Jacobian at the equilibrium $E_{0}$:

[TABLE]

where

[TABLE]

and

[TABLE]

Then the characteristic equation of the Jacobian matrix $J_{E_{0}}$ is

[TABLE]

where the variable $x$ denotes eigenvalues of the matrix $J_{E_{0}}$. Obviously, when $\mathcal{R}_{0}<1$, all the eigenvalues are negative and the DFE point $E_{0}$ is locally asymptotically stable. Otherwise, when $\mathcal{R}_{0}>1$, there is a positive eigenvalue, suggesting that the DFE point $E_{0}$ is unstable. This concludes the proof. ∎

Furthermore, in the following theorem we provide a sufficient condition to guarantee the global asymptotically stability of DFE.

Theorem 4**.**

If $\mathcal{R}_{0}<1$ and $\lambda\leq\varphi$, then the DFE point $E_{0}$ of model (5) is globally asymptotically stable.

Proof 3**.**

Suppose $\lambda\leq\varphi$, then it follows from the first equation of model (5) that for any $r=r_{1},r_{2},\dots,r_{m}$, we have

[TABLE]

Hence, $s(r,t)\leq\frac{\varphi}{\varphi+\theta(r)}$ for any $r=r_{1},r_{2},\dots,r_{m}$ and $t\geq 0$. Define the Lyapunov function as

[TABLE]

Then, the derivative of $L(t)$ along the system (5) is given as

[TABLE]

Note that $\mathcal{R}_{0}<1$ is equivalent to $\beta\rho\varphi\pi{\sum_{j=1}^{m}P(r_{j}){\frac{r_{j}^{2}}{\varphi+\theta(r_{j})}}}-\lambda<0$. Therefore, if $\lambda\leq\varphi$ and $\mathcal{R}_{0}<1$, then $L^{\prime}(t)\leq 0$. Furthermore, $L^{\prime}(t)=0$ if and only if $i(r,t)=0$ for any $r=r_{1},r_{2},\dots,r_{m}$. By LaSalle’s invariance principle [74], we conclude that if $\mathcal{R}_{0}<1$ and $\lambda\leq\varphi$, then the DFE point $E_{0}$ is globally asymptotically stable. This completes the proof. ∎

4.2 Stability of EE

Theorem 5**.**

If $\mathcal{R}_{0}>1$, then the system (5) admits a unique EE point $E^{*}$ defined by Eq. (12) which is locally asymptotically stable in $\Gamma$.

Proof 4**.**

Define $\hat{s}(r,t)=s(r,t)-s_{r}^{*}$ and $\hat{i}(r,t)=i(r,t)-i_{r}^{*}$ for each $r=r_{1},r_{2},\dots,r_{m}$. We consider the following linearized dynamics of system (5) at $E^{*}$.

[TABLE]

where $\hat{I}(t)=\sum\limits_{j=1}^{m}P(r_{j})\hat{i}(r_{j},t)$ and $I^{*}=\sum\limits_{j=1}^{m}P(r_{j})i_{r_{j}}^{*}$.

Suppose $\xi$ is an arbitrary eigenvalue of the coefficient matrix of the linearized system (13). Then the proof can be done if we verify that $\xi$ has negative real part, i.e., $Re{(\xi)}<0$. Looking for exponential solutions of the linear equations (13), we set $\hat{s}(r,t)=\hat{s}_{0}(r)e^{\xi t}$ and $\hat{i}(r,t)=\hat{i}_{0}(r)e^{\xi t}$ for any $r=r_{1},r_{2},\dots,r_{m}$. Substituting in the linearized system and canceling $e^{\xi t}$, we obtain

[TABLE]

Thus

[TABLE]

By denoting $Y(\hat{i}_{0})=\sum\limits_{j=1}^{m}P(r_{j})\hat{i}_{0}(r_{j})=\hat{I}(t)e^{-\xi t}$, we get

[TABLE]

Multiplying Eq. (16) by $P(r)e^{\xi t}$ and summarizing over all $r=r_{1},r_{2},\dots,r_{m}$ gives rise to

[TABLE]

If $\hat{I}(t)=0$, then $Y(\hat{i}_{0})=0$ and hence $\hat{i}_{0}(r)=0$ by Eq. (16). In this case, we obtain from the first equation of (14) that

[TABLE]

If $\hat{I}(t)\neq 0$, then according to Eq. (17) we have

[TABLE]

Letting the right-hand side of the second equation of system (5) be zero, we obtain a stationary solution

[TABLE]

By Eq. (9) we obtain

[TABLE]

Note that $I^{*}>0$ as long as $\mathcal{R}_{0}>1$, therefore

[TABLE]

Next, we verify that $Re(\xi)<0$ by the way of contradiction.

(i) if $\xi=0$, then it follows from Eq. (18) that

[TABLE]

Since $I^{*}>0$, Eq. (20) indicates that

[TABLE]

which is a contradiction. Hence, $\xi\neq 0$.

(ii) if $Re(\xi)>0$, then $Re(\xi+\frac{\xi+\varphi}{\xi+\theta(r_{j})+\varphi}\beta\rho\pi r_{j}^{2}I^{*})>0$ $(j=1,2,\dots,m)$. Let

[TABLE]

where $\mathbf{i}$ is the imaginary unit. Consequently, we have

[TABLE]

Therefore,

[TABLE]

This implies

[TABLE]

which gives a contradiction to Eq. (18). Hence, $Re(\xi)\leq 0$.

(iii) If $Re(\xi)=0$, then the imaginary part $Im(\xi)\neq 0$ according to (i). In this case, we also find that $Re(\xi+\frac{\xi+\varphi}{\xi+\theta(r_{j})+\varphi}\beta\rho\pi r_{j}^{2}I^{*})>0$ $(j=1,2,\dots,m)$. This results in a contradiction with Eq. (18) again in an analogous fashion. Thus, $Re(\xi)<0$.

Based on the above discussion, we conclude that if $\mathcal{R}_{0}>1$, then all the eigenvalues have negative real part, thus the EE point $E^{*}$ of model (5) is locally asymptotically stable. This completes the proof. ∎

5 Persistence of the disease

As pointed out in the previous section, there is an EE as long as $\mathcal{R}_{0}>1$. In this section we additionally present the following theorem on the persistence of the disease in the case of $\mathcal{R}_{0}>1$.

Theorem 6**.**

If $\mathcal{R}_{0}>1$, then the system (5) is persistent, that is, there exists $\varepsilon>0$ such that

[TABLE]

Proof 5**.**

We will use the conclusion given by Thieme (see theorem 4.6 in [75]) to prove the above proposition. Starting with the positively invariant set $\Gamma$ given by (6), we define two sets

[TABLE]

According to the proof of Lemma 1, if $s(r,0)\geq 0,v(r,0)\geq 0,I(0)=\sum\limits_{j=1}^{m}P(r_{j})i(r_{j},0)>0$, then $s(r,t)\geq 0,v(r,t)\geq 0,I(t)=\sum\limits_{j=1}^{m}P(r_{j})i(r_{j},t)>0$ for any $t>0$. Note that

[TABLE]

and $I(0)=\sum\limits_{j=1}^{m}P(r_{j})i(r_{j},0)>0$, it follows from the comparison theorem [76] that

[TABLE]

Hence, the set $\Gamma^{*}$ is also positively invariant. Moreover, there exists a compact set $B$ in which all solutions to (5) starting from $\Gamma$ will enter and remain permanently. It can be easily validated that the set $B$ satisfies the compactness conditions ($C_{4.2}$) proposed by Thieme [75]. Denote

[TABLE]

and

[TABLE]

where \omega\Big{(}s(r,0),i(r,0)\Big{)} is the $\omega$-limit set of the solutions of system (5) initiated from $(s(r,0),i(r,0))$. Confining the system (5) to $W_{\partial}$ yields

[TABLE]

It is easy to manifest that the system (21) has a unique equilibrium $E_{0}$ as denoted by Eq. (7), which is globally asymptotically stable. Thus, $\Omega_{\Gamma}=\{E_{0}\}$. In addition, $E_{0}$ is an acyclic isolated covering of $\Omega_{\Gamma}$ since there is no solution in $W_{\partial}$ that links $E_{0}$ to itself. Next, we will prove that $\{E_{0}\}$ is a weak repeller for $\Gamma^{*}$, namely, any solution $(s(r,t),i(r,t))$ with initial value in $\Gamma^{*}$ satisfies

[TABLE]

Here, the distance $dist(x,Y)$ of a point $x\in X$ from a subset $Y$ of $X$ is defined by

[TABLE]

where $d$ is a metric of the metric space $X$. According to the proof of lemma 3.5 by Leenheer and Smith [77], we only need to verify $M^{s}(E_{0})\cap\Gamma^{*}=\emptyset$, where $M^{s}(E_{0})$ is the stable manifold of $E_{0}$. We prove it by the way of contradiction. Suppose $M^{s}(E_{0})\cap\Gamma^{*}\neq\emptyset$, then there exists a solution $(s(r,t),i(r,t))\in\Gamma^{*}$ such that

[TABLE]

It is clear form (22) that for any given $\eta>0$, there exists $T_{0}>0$ such that

[TABLE]

Given that $\mathcal{R}_{0}=\frac{\beta\rho\pi}{\lambda}\sum\limits_{j=1}^{m}P(r_{j})r_{j}^{2}{\frac{\varphi}{\varphi+\theta(r_{j})}}>1$, there exists a positive constant $\eta>0$ such that

[TABLE]

Let

[TABLE]

thus $C>0$ as long as $\mathcal{R}_{0}>1$. Define the following Lyapunov function

[TABLE]

then the derivative of $V$ along the solution $(s(r,t),i(r,t))$ of system (5) is given by

[TABLE]

It follows that for all $t\geq T_{0}$,

[TABLE]

That is, \frac{dV(t)}{dt}\Big{|}_{\eqref{func4}}\geq CV(t) for all $t\geq T_{0}$, which implies $\lim\limits_{t\to\infty}V(t)=\infty$. This is a contradiction to the boundedness of $V(t)$. Therefore, $\{E_{0}\}$ is a weak repeller for $\Gamma^{*}$. Based on the result by Thieme (theorem 4.6 in [75]), we conclude that if $\mathcal{R}_{0}>1$ then the system (5) is persistent. ∎

6 Simulation results and discussion

In this section, we provide extensive stochastic simulations to support the theoretical results of our model. In the end we also display some simulation results on the topological structure of the underlying contact network.

6.1 Spreading dynamics of the epidemic model

In our simulations, the length of the side of the square space is set to be $D=30$ and the velocity $v$ of each individual is fixed to be $v=0.1$. We start the simulation with $1\%$ of individuals being infected seeds while the others all susceptible. That is, at $t=0$ the number of infected, susceptible and vaccinated individuals is $N_{I}(0)=9$, $N_{S}(0)=891$ and $N_{V}(0)=0$, respectively. At the beginning, all the individuals are randomly distributed within the planar space, where at each time step the individuals move depending on the random jump probability $p_{\text{j}ump}$. As long as an infected individual runs into an effective interaction radius of a susceptible one, an infection takes place at the transmission rate $\beta$. The simulation is ended as the number (or fraction) of each class of individuals reaches a relatively steady level (with negligible fluctuations). To investigate the impacts of heterogeneity in the interaction radius among individuals, we consider different distributions of interaction radius in the model, including the Poisson, exponential, power-law and in some cases the Kronecker delta.

Figure 4 exhibits a comparison of the value of the basic reproduction number $\mathcal{R}_{0}$ as a function of the vaccination-strength-related parameter $\alpha$ between four different distributions $P(r)$ of interaction radius of individuals. In the case of the Kronecker Delta distribution $P(r)=\delta_{r3}$ (where $\delta_{r3}=1$ if $r=3$ and $\delta_{r3}=0$ otherwise), all the individuals have an identical interaction radius, thus the vaccination rate $\theta(r)$ turns to be $\theta(r)=\theta(3)=\theta_{0}$ which is a constant. This indicates the value of $\mathcal{R}_{0}$ is a constant, in particular by Eq. (11), $\mathcal{R}_{0}=3^{2}\beta\varphi\rho\pi/[\lambda(\varphi+\theta_{0})]\simeq 0.135$. In the case of Poisson, exponential and power-law distributions, the value of $\mathcal{R}_{0}$ remains almost unchanged for $\alpha<-5$. This is because the smaller the value of parameter $\alpha$, the larger the vaccination probability for individuals with small radius. Therefore, as $\alpha$ is small enough, the individuals to be vaccinated are only those with the smallest radius since their vaccination probability is dominantly large. When $-5<\alpha<2$, the value of $\mathcal{R}_{0}$ decreases with $\alpha$ drastically and in general the value of $\mathcal{R}_{0}$ in the power-law radius distribution is larger than the value of $\mathcal{R}_{0}$ in the exponential radius distribution, which is in turn greater than that in the Poisson radius distribution. When $\alpha>2$, the value of $\mathcal{R}_{0}$ first increases relatively fast and then grows gradually slowly with large $\alpha$. Again, the values of $\mathcal{R}_{0}$ in power-law and exponential radius distributions are greater than that in the Poisson and Delta radius distributions. This result suggests that it is easier for the disease to break out in the population with more heterogeneous distribution of interaction radius due to a larger $\mathcal{R}_{0}$. It is interesting to notice that for all of the power-law, exponential and Poisson distributions of interaction radius, the value of $\mathcal{R}_{0}$ reaches a minimum at $\alpha=2$, as illustrated by the grey vertical line in Fig. 4. This may motivate an optimal vaccination intervention for disease prevention irrespective of the distribution of individuals’ effective interaction radius.

Figure 5 depicts the time series of the density of each class of individuals in the population for different values of $p_{\text{j}ump}$ given other parameters. It is clear from Fig. 5(b) that the density of infected individuals grows with $p_{\text{j}ump}$, while it is shown from Figs. 5(a, c) that both the densities of susceptible individuals and vaccinated ones decrease with $p_{\text{j}ump}$. This means that for larger $p_{\text{j}ump}$ the infection is more severe [35, 36, 37, 50]. Moreover, the simulation results in the case with $p_{\text{j}ump}=1$ are in perfect agreement with the theoretical predictions based on the HM assumption, which confirms the statement ahead of model (2). As done in Ref. [50], in what follows we only investigate the case of $p_{\text{j}ump}=1$ since in other cases our model (2) based on HM assumption deviates obviously from the simulation results, as demonstrated in Fig. 5.

Figure 6 plots the temporal evolution of the density of susceptible, infected and vaccinated individuals in the population for different distributions of effective interaction radius. Both the results under the condition of $\mathcal{R}_{0}<1$ and in the case of $\mathcal{R}_{0}>1$ have been calculated. When $\mathcal{R}_{0}<1$, the fraction of infected population eventually decays to zero irrespective of the radius distribution. Otherwise, when $\mathcal{R}_{0}>1$, as shown in Figs. 6(b,d,f) the density of infected individuals raises swiftly in the early stage and finally enters a stationary state. It is demonstrated that all the simulation results in the case of $p_{\text{j}ump}=1$ agree well with the theoretical predictions by our model based on the HM assumption. To further explore the effects of the heterogeneity in the interaction radius on the epidemic spread, we give a clear comparison of the density of infected individuals between the Poisson, exponential and power-law radius distributions in Fig. 7, where the data is extracted from Figs. 6(b,d,f). It is shown from Fig. 7 that the final epidemic prevalence in the power-law interaction radius distribution is at the lowest level although the basic reproduction number $\mathcal{R}_{0}=4.655$ in the power-law case is the largest among the three distributions. On the contrary, the final epidemic prevalence (i.e. the final fraction of infected individuals) in the case of Poisson radius distribution is the highest, albeit with the smallest basic reproduction number $\mathcal{R}_{0}=3.5017$. This interesting result reveals that the disease will break out readily in the population with a heterogeneous interaction radius distribution; however, resulting in a relatively small epidemic prevalence.

Figure 8 displays the dependence of the final densities of individuals in each class on the tunable parameter $\alpha$ that reflects the vaccination strength relevant to the individual’s effective interaction radius (1). In the left column panels the results are sorted out according to different distributions of interaction radius: (a) Poisson, (b) exponential and (c) power-law with the same average radius $\langle r\rangle=3$. Conversely, in the right column panels the same results are sorted out according to different disease status: (d) susceptible, (e) infected and (f) vaccinated. It is shown that the final density of susceptible (vaccinated) individuals reaches the minimum (maximum) at about $\alpha=0$, without respect to the radius distribution (see Figs. 8(d,f)). In the case of $\alpha<0$, the vaccination favors susceptible individuals with small interaction radius. Moreover, the smaller the parameter $\alpha$, the larger the vaccination probability for susceptible individuals with smaller radius. Consequently, the final density of susceptible (vaccinated) individuals increases (decreases) as $\alpha$ decays since the vaccination covers only the susceptible individuals whose interaction radius is relatively small. As $\alpha$ becomes small enough, the vaccination only covers the susceptible individuals with the smallest interaction radius. Therefore the final densities of susceptibles and vaccinateds remain constant when $\alpha$ is extremely small. In the case of $\alpha>0$, the vaccination favors susceptible individuals with large interaction radius. The larger the value of $\alpha$, the higher vaccination probability for susceptible individuals with larger interaction radius. This implies that the final density of susceptible (vaccinated) individuals increases (decreases) with the increase of $\alpha$ since the vaccination takes place only for the susceptible individuals whose interaction radius is relatively large. However, when $\alpha$ is big enough, the vaccination only covers the susceptible individuals with the biggest interaction radius. Hence the final density of susceptible (vaccinated) individuals remains almost constant as $\alpha$ is extremely large. In general, as the value of $\alpha$ grows gradually, the final infected density first keeps almost unchanged, then drops to zero and remains in the disease-free state for a range of $\alpha$ with $\mathcal{R}_{0}<1$, and then rises fast and finally grows very slowly. The differences of the final density of each class of individuals between different radius distributions have also been demonstrated in Figs. 8(d,e,f). In particular, as the value of $\alpha$ is large enough or small enough, the final density of infected individuals for the power-law radius distribution is lower than the infected density in the exponential radius distribution, which in turn is lower than that in the Poisson radius distribution (see Fig. 8(e)). Otherwise, the interval of $\alpha$ that satisfies $\mathcal{R}_{0}<1$ for the Poisson radius distribution is wider than that in the exponential and power-law radius distributions. All the above results illustrate a good agreement between stochastic simulations and theoretical calculations, additionally the extreme points at about $\alpha=0$ may suggest an optimal vaccination intervention for disease control, regardless of the radius distribution.

We also extract the final epidemic prevalence as a function of the density $\rho$ of moving individuals. Results are reported in Fig. 9 where a comparison is exhibited among different distributions of interaction radius with the same expectation $\langle r\rangle=3$ and the same parameters of infection. It is observed that given all other parameters there exists a critical value of density $\rho_{\rm c}$, above which the disease breaks out, otherwise it dies out. In fact it can be deduced from Eq. (11) that the threshold condition $\rho=\rho_{\rm c}$ is equivalent to $\mathcal{R}_{0}=1$, and that $\rho>\rho_{\rm c}$ $(\rho<\rho_{\rm c})$ is equivalent to $\mathcal{R}_{0}>1$ $(\mathcal{R}_{0}<1)$. On the other hand, when $\rho>\rho_{\rm c}$ the epidemic prevalence versus $\rho$ is monotonic, as intuition suggests that individuals in a population with higher density are more connected to each other, leading to a higher level of infection. In addition, it is shown that the value of $\rho_{\rm c}$ under the power-law and exponential radius distributions is smaller than that with Poisson radius distribution. As $\rho$ is far larger than $\rho_{\rm c}$ the epidemic prevalence in the population with Poisson radius distribution is larger than that with exponential radius distribution, which is in turn higher than that in the population with power-law radius distribution. This result restates the conclusion drawn from Fig. 7 that the disease inclines to spread in the population with more heterogeneous interaction radius; however, resulting in a relatively lower epidemic prevalence.

6.2 Structure of the directed contact network

Figure 10 presents both the in-degree distribution $P(k_{\rm in})$ and the out-degree distribution $P(k_{\rm out})$ of the resulting directed contact network in the steady state for different interaction radius distributions of individuals. As shown in Fig. 10(a), for moving individuals with power-law (marked by pink diamonds) or exponential (marked by blue triangles) radius distributions, the in-degrees of the directed contact network follow a multimodal distribution with a power-law decaying trend as shown in the log-log plot. In the case of Kronecker Delta radius distribution (marked by red squares), the in-degrees of the directed contact network follows a Poisson distribution with an expectation of $\langle k_{\rm in}\rangle=28.5$ (see the left inset in Fig. 10(a)), while in the case of Poisson radius distribution (marked by black circles), the in-degree of the directed contact network obeys a multimodal distribution with an exponential decaying trend as shown in the semilog plot (see the right inset in Fig. 10(a)). On the contrary, Fig. 10(b) indicates that the out-degrees of the directed contact network follow a Poisson distribution, in which the expectation $\langle k_{\rm out}\rangle$ is relevant to the radius distribution. As illustrated by the Poisson fitting lines in Fig. 10(b), the average out-degree $\langle k_{\rm out}\rangle$ of the directed network is $28.5,35.5,42.6,47.4$ for Delta, Poisson, exponential and power-law distributions of interaction radius, respectively. The more heterogeneous the interaction radius distribution, the larger the average out-degree of the contact network.

In Fig. 11 we provide the average in-degree $k_{\rm in}(r)$ and out-degree $k_{\rm in}(r)$ of individuals that have radius $r$ as a function of the radius $r$. It is clear from Fig. 11(a) that $k_{\rm in}(r)$ is proportional to the area $\pi r^{2}$ of the effective contact space (circle) of individuals with radius $r$. This relation is not surprising and it is reasonable since the in-degree of an individual is defined as the total number of other individuals that fall within the circle with respect to the radius of the individual considered. Therefore, the larger the area of the contact circle, the larger the in-degree of the individuals. In this case, note that the density is $\rho=N/D^{2}=1$, we have exactly $k_{\rm in}(r)=\pi r^{2}$, as indicated by the black solid line in Fig. 11(a). Moreover, it is shown that this result is independent on the radius distribution. In a sharp contrast, it is observed from Fig. 11(b) that the average out-degree $k_{\rm out}(r)$ of individuals with radius $r$ is completely uncorrelated with the radius $r$. However, the value of the average out-degree is strongly related to the radius distribution. As demonstrated in Fig. 11(b), we have $k_{\rm out}(r)=35.5$, $42.6$ and $47.4$ for Poisson, exponential and power-law distributions of interaction radius, respectively.

Figure 12 outlines the correlation between in-degree $k_{\rm in}$ and out-degree $k_{\rm out}$ for different cases of radius distributions. This correlation is determined by calculating the average in-degree of individuals whose out-degree is $k_{\rm out}$. In the cases of Poisson, exponential and power-law radius distributions, the correlation between $k_{\rm in}$ and $k_{\rm out}$ takes the unimodal form. As the out-degree increases, the in-degree first increases towards a summit and then decreases gradually to zero. In the case of the Delta radius distribution, we observe a linear correlation with $k_{\rm in}$ being equal to $k_{\rm out}$.

In order to deep understand how the vaccination behavior affects the network structure and the transmission process, in Fig. 13 the average in-degree and average out-degree of individuals in each class are plotted against the vaccination-strength-related parameter $\alpha$ for different distributions of interaction radius. It is shown from the Figs. 13(a,c,e) that the average in-degrees $\langle k_{\rm in}^{S}\rangle$, $\langle k_{\rm in}^{I}\rangle$ and $\langle k_{\rm in}^{V}\rangle$ of susceptible, infected and vaccinated individuals reached their respective extreme values at $\alpha=2$, irrespective of the radius distribution. A direct comparison with Fig. 4 suggests that the average in-degree of individuals is implicitly relevant to the behavior of the basic reproduction number $\mathcal{R}_{0}$ which arrives at the minimal value at $\alpha=2$. On the other hand, Figs. 13(b,d,f) show that the behavior of the average out-degree of individuals in each class closely resembles the final density of each class of individuals (see Fig. 8) in that they reach their extreme values at about $\alpha=0$ without regard to the radius distribution. This similarity indicates that the average out-degree of individuals is potentially responsible for the behavior of the final density with regard to $\alpha$.

7 Conclusions

In summary, we have established an SIS epidemic model with vaccination on a dynamic network of mobile individuals with spatial constraints where all individuals have heterogeneous interaction radii (or contact circles) in a two-dimensional space. In the model, we consider that the vaccination of each susceptible individual depends on the individual’s interaction radius and assume that all individuals are random walkers who are also allowed to perform a long-distance jump with a probability $p_{\rm jump}$. We derive a homogeneous mixing model with a set of ordinary differential equations in the special case of $p_{\rm jump}=1$. We have obtained the basic production number $\mathcal{R}_{0}$ and studied the dynamical behavior of the model. We argue that the disease-free equilibrium is locally asymptotically stable as $\mathcal{R}_{0}<1$; otherwise, if $\mathcal{R}_{0}>1$ then the disease-free equilibrium is unstable and there is a unique endemic equilibrium which is locally asymptotically stable and persistent. Extensive computational simulations have been carried out to further explore the dynamical behavior of the epidemic spreading as well as the topological properties of the underlying contact network. Our results have demonstrated a good agreement between theory and simulations with regard to the disease transmission. It is found that different distributions of individuals’ interaction radius have strong impacts on the basic reproduction number and the final densities of individuals of each class. Generally speaking, the heterogeneity of interaction radius among individuals will facilitate the disease transmission while resulting in a relatively low epidemic prevalence. Moreover, based on the dependencies of $\mathcal{R}_{0}$ and the final epidemic prevalence on the vaccination-strength-related parameter $\alpha$, we argue that an optimal vaccination intervention is achievable for disease prevention and control. Furthermore, it is shown that the in-degree distribution of the resulting network in the cases of power-law and exponential radius distributions follows a multimodal distribution with a power-law decaying trend, whereas in the case of Poisson radius distribution, the in-degrees of the contact network follows a multimodal distribution with an exponential decaying trend. Conversely, in the case of Delta radius distribution, the in-degree distribution is Poisson. On the other hand, the out-degrees of the contact network follows a Poisson distribution with different expectations for different radius distributions. The correlation between in-degree and out-degree of the contact network takes the form of a unimodal function. It is worth noticing that the average in-degree and the average out-degree of individuals in each class have the same qualitative behavior with respect to the parameter $\alpha$. This resemblance provides another perspective to understand the influence of the vaccination intervention on the epidemic transmission as it has great impacts on the average in-degree and out-degree of individuals. This work provides an in-depth analysis for the dynamical behavior of the epidemic model and sheds new light on potentially optimal vaccination interventions for epidemic spreading in the population of moving individuals with spatial limitations.

Acknowledgments

This research project was jointly supported by the National Natural Science Foundation of China (Grant Nos. 11601294, 61873154, 11501340, 11331009), Shanxi Province Science Foundation for Youths (Grant No. 201601D021012) and Shanxi Scholarship Council of China (Grant No. 2016-011).