D-optimal design of the additive mixture model with multi-response

1.   Introduction

Population migration is a common phenomenon. With the migration of population, infectious diseases can easily spread from one area to another, so it is meaningful to consider population migration when studying the spread of infectious diseases ^{[1,2,3,4,5,6,7]}.

Wang and Mulone ^[2] established an SIS infectious disease model with standard incidence based on two patches. It is proved that the basic reproduction number is the threshold of the uniform persistence and disappearance of the disease. The dispersal rate of the population will make the infectious disease persist or disappear in all patches. There will be no the phenomena that infectious diseases persists in one patch but disappears in the other.

Sun et al. ^[3] put forward an SIS epidemic model with media effect in a two patches setting. Under the assumption that the migration matrix is irreducible, it is proved that if the basic reproduction number is greater than 1 then the system persists and solutions converge to an endemic equilibrium and that if the basic reproduction number is less than 1 then solutions tend to an equilibrium without disease.

Gao et al. ^[5,6] studied an SIS multi-patch model with variable transmission coefficients. Their results show that the basic reproduction number $R_0$ is a threshold parameter of the disease dynamics.

All the patch models referenced above assume that the migration matrix is irreducible. The studies in which the migration matrix is reducible are few. Therefore, based on the case of two patches, we consider that the individuals can only migrate from one patch to the other. In this case, the migration matrix is reducible. It can characterize the phenomenon that individuals migrate in one direction between two regions, such as, from the rural patch to the urban one ^[8] and from a small community hospital to a large teaching hospital ^[4].

It is well known that the incidence rate plays an important role in the modeling of infectious disease. Considering the saturation phenomenon for numerous infected individuals, Capasso and Serio ^[9] first introduce a nonlinear bounded function $g(I)$ to form the interaction term $g(I) S$ in 1978. It can characterize the behavioral changes of individuals, such as wearing masks or reducing their social activities and direct contact with others with the increase of infectious individuals. After that, the saturation incidence rate has attracted much attention and various nonlinear types of incidence rate are employed. The most commonly used types are Holling type Ⅱ $\frac{\lambda SI}{1+\alpha I}$ ^[10,11,12] and $\frac{\beta SI}{1+\alpha S}$ ^[13], Monod-Haldane type $\frac{\lambda SI}{1+\alpha I^2}$ ^[14], Beddington-DeAngelis type $\frac{\lambda SI}{1+\alpha S+\beta I}$ ^[15,16,17] and Crowley-Martin type $\frac{\lambda SI}{(1+\alpha S)(1+\beta I)}$ ^[18,19].

In this paper, we consider infectious disease transmission models with saturation incidence rate. The rest of this paper is organized as follows: In Section 2, we establish a two-patch SIS model with saturating contact rate and one-directing population dispersal. We discuss the existence of disease-free equilibrium, boundary equilibrium and endemic equilibrium and prove the global asymptotic stability of the equilibriums in Section 3. In Section 4, we perform simulations to illustrate the results and analyze the effect of the contact rate and population migration on epidemic transmission. Finally, we discuss in Section 5.

2.   Model formulation

In the two patches, the population is divided into two states: susceptible and infective. Thus we can establish a two-patch SIS model with saturating contact rate and one-directing population dispersal

where $S_i$ is the number of susceptible population in patch $i$ $(i = 1, 2)$, $I_i$ is the number of infective population in patch $i$ $(i = 1, 2)$, $A_i$ is the recruitment into patch $i$ $(i = 1, 2)$, $d_i$ is the natural mortality rate, $\gamma_i$ is the recovery rate of an infective individual in patch $i$ $(i = 1, 2)$, $m$ is the migration rate form patch 1 to patch 2. Since the individuals can migrate from the first patch to the second, patch 1 is the source patch and patch 2 is the sink patch. The initial conditions is

Denote the population in patch $i$ by $N_i$. Then $N_i = S_i+I_i$. From system (2.1), the differential equations governing the evolution of $N_1$ and $N_2$ are

Obviously, system (2.3) has a unique equilibrium $(N_1^*, N_2^*) = (\frac{A_1}{d_1+m}, \frac{A_2}{d_2}+\frac{mA_1}{d_2(d_1+m)})$ which is globally asymptotically stable for (2.3). So (2.1) is equivalent the following system

Because $\lim\limits_{t\rightarrow \infty}N_i(t)\rightarrow N_i^*$ $(i = 1, 2)$, system (2.4) leads to the following limit system

3.   Mathematical analysis for system (2.5)

3.1.   The invariants and equilibriums

Let $\Omega = \{(I_1, I_2)|0\leq I_1\leq N_1^*, 0\leq I_2\leq N_2^*\}$. Then $\Omega$ is invariant region for system (2.5).

Define the basic reproduction number in the two patches respectively by $R_0^1 = \frac{\beta_1A_1}{(d_1+\gamma_1+m)^2} = \frac{\beta_1}{d_1+\gamma_1+m}N_1^*$, $R_0^2 = \frac{\beta_2}{(d_2+\gamma_2)}N_2^*$. The basic reproduction number $R_0^1$ gives the expected secondary infections in the source patch produced by a primary infected individual in the source patch when the population is supposed to be in the disease-free equilibrium. The basic reproduction number $R_0^2$ gives the expected secondary infections in the sink patch produced by a primary infected individual in the sink patch when the population is supposed to be in the disease-free equilibrium. Then we have the following theorem.

Theorem 3.1. For the system (2.5), we have

(i) The disease-free equilibrium $E_0: = (0, 0)$ always exists;

(ii) The boundary equilibrium $E_1: = (0, \frac{\beta_2 N_2^*-d_2-\gamma_2}{(d_2+\gamma_2)\alpha_2+\beta_2})$ exists if $R_0^2 > 1$;

(iii) There is a unique epidemic equilibrium $E_*$ if $R_0^1 > 1$.

Proof. (ⅰ) can be easily proved.

Let

From Eq (3.1), we can have $I_1 = 0$ always satisfies Eq (3.1). When $I_1 = 0$, from Eq (3.2), we have

If $R_0^2 > 1$, then $I_2 = \frac{\beta_2N_2^*-d_2-\gamma_2}{(d_2+\gamma_2)\alpha_2+\beta_2} > 0$. So The boundary equilibrium $E_1: = (0, \frac{\beta_2N_2^*-d_2-\gamma_2}{(d_2+\gamma_2)\alpha_2+\beta_2})$ exists if $R_0^2 > 1$. The conclusion (ⅱ) is proved.

If $R_0^1 > 1$, Eq (3.1) has a positive solution $I_1^* = \frac{\beta_1 N_1^*-(d_1+m+\gamma_1)}{(d_1+m+\gamma_1)\alpha_1+\beta_1}$. Solve Eq (3.2), we have

Substituting $I_1^*$ into Eq (3.3), we have

Since $I_2^*\geq 0$ is meaning only, we take

So if $R_0^1 > 1$, there is a unique epidemic equilibrium $E_* = (I_1^*, I_2^*)$, where $I_1^* = \frac{\beta_1 N_1^*-(d_1+m+\gamma_1)}{(d_1+m+\gamma_1)\alpha_1+\beta_1}$ and $I_2^* = \frac{(\beta_2 N_2^*-d_2-\gamma_2+m\alpha_2 I_1^*)+\sqrt{(\beta_2 N_2^*-d_2-\gamma_2+m\alpha_2 I_1^*)^2+4[(d_2+\gamma_2)\alpha_2+\beta_2]mI_1^*}}{2[(d_2+\gamma_2)\alpha_2+\beta_2]}.$ The conclusion (ⅲ) is proved.

This completes the proof of the theorem.

From the above analysis, we have the following theorem.

Theorem 3.2. For the system (2.5), we have

(i) If $R_0^1 < 1$ and $R_0^2 < 1$, there is the disease-free equilibrium $E_0$ only;

(ii) If $R_0^1 < 1$ and $R_0^2 > 1$, there are the disease-free equilibrium $E_0$ and the boundary equilibrium $E_1$;

(iii) If $R_0^1 > 1$ and $R_0^2 < 1$, there are the disease-free equilibrium $E_0$ and the epidemic equilibrium $E_*$;

(iv) If $R_0^1 > 1$ and $R_0^2 > 1$, there are the disease-free equilibrium $E_0$, the boundary equilibrium $E_1$ and the epidemic equilibrium $E_*$.

Remark 3.1. Define the basic reproduction number $R_0$ of the system (2.5) by the spectral radius of the next generation matrix ^[20], we have

where $\rho(A)$ denotes the spectral radius of a matrix $A$. So from the above analysis, we know that $R_0 = \max\{R_0^1, R_0^2\}$.

3.2.   The stability of the disease-free equilibrium

The next, we shall discuss the local stability of the disease-free equilibrium firstly. Then we discuss the global asymptotical stability.

Theorem 3.3. For the system (2.5), we have

(i) If $R_0^1 < 1$ and $R_0^2 < 1$, the disease-free equilibrium $E_{0}$ is locally asymptotically stable;

(ii) If $R_0^1 > 1$ or $R_0^2 > 1$, the disease-free equilibrium $E_{0}$ is unstable.

Proof. The linearized system of (2.5) at the equilibrium $E_{0}$ is

The associated characteristic equation of the linearized system of (3.4) at the equilibrium $E_0$ is

It is easy to see that the two eigenvalues of characteristic Eq (3.5) are

and

So, when $R_0^1 < 1$ and $R_0^2 < 1$, the disease-free equilibrium $E^{0}$ is locally asymptotically stable; However, if $R_0^1 > 1$ or $R_0^2 > 1$, the disease-free equilibrium $E^{0}$ is unstable.

Remark 3.2. From Theorem 3.3, we know that for the system (2.5), if $R_0 < 1$ the disease-free equilibrium $E_{0}$ is locally asymptotically stable; if $R_0 > 1$, the disease-free equilibrium $E_{0}$ is unstable.

Theorem 3.4. For the system (2.5), if $R_0^1 < 1$ and $R_0^2 < 1$, the disease-free equilibrium $E_{0}$ is globally asymptotically stable.

Proof. Since $\frac{I_i}{1+\alpha_iI_i}\leq I_i$ for $i = 1, 2$, from system (2.5), we can obtain that

Define an auxiliary linear system using the right hand side of (3.6) as follows

It can be rewritten as

if $R_0^1 < 1$ and $R_0^2 < 1$, we can solve (3.7) and know that $\lim\limits_{t\rightarrow \infty}I_1(t) = 0$ and $\lim\limits_{t\rightarrow \infty}I_2(t) = 0$. By the comparison principle ^[21], we can conclude that when $R_0^1 < 1$ and $R_0^2 < 1$, all non-negative solutions of (2.5) satisfy $\lim\limits_{t\rightarrow \infty}I_i(t) = 0$ for $i = 1, 2$. So the disease-free equilibrium $E_{0}$ is globally asymptotically stable.

3.3.   The stability of the boundary equilibrium

In this subsection, we will discuss the local stability of the boundary equilibrium firstly. Then discuss the global asymptotical stability.

Theorem 3.5. For the system (2.5), if $R_0^1 < 1$ and $R_0^2 > 1$, the boundary equilibrium $E_{1}$ is globally asymptotically stable.

Proof. The Jacobian matrix at the boundary equilibrium $E_{1}$ of system (2.5) is

The two eigenvalues of the Jacobian matrix are

and

So, when $R_0^1 < 1$ and $R_0^2 > 1$, $\lambda_1 < 0$ and $\lambda_2 < 0$. That is the boundary equilibrium $E_{1}$ is locally asymptotically stable.

For every $(I_1(0), I_2(0)) \in \Omega$, assume the solution of the system (2.5) with initial value $(I_1(0), I_2(0))$ is $(I_1(t), I_2(t))$. Since

if $R_0^1 < 1$, $\dfrac{dI_1(t)}{dt} < 0$, then $I_1(t)$ is positive and decreasing and $\lim\limits_{t\rightarrow \infty}I_1(t) = 0$. So for sufficiently small positive number $\epsilon_1$, there exists a T, such that $I_1(T) = \epsilon_1$ and $I_1(t) < \epsilon_1$ when $t > T.$

The following, we prove that for any $\epsilon > 0$, there exists a $T^* > T$ such that $|I_2(T^*)-\frac{\beta_2N_2^*-d_2-\gamma_2}{\beta_2+(d_2+\gamma_2)\alpha_2}| < \epsilon$. And because $E_{1} = (0, \frac{\beta_2N_2^*-d_2-\gamma_2}{\beta_2+(d_2+\gamma_2)\alpha_2})$ is locally asymptotically stable, we have $E_{1}$ is globally asymptotically stable.

Since

if $I_2(T) < \frac{\beta_2N_2^*-d_2-\gamma_2}{\beta_2+(d_2+\gamma_2)\alpha_2}$, then $\dfrac{dI_2(t)}{dt} > 0$ for $t > T.$ So $I_2(t)$ is increasing and there exists $T^*_1$ such that $|I_2(T^*_1)-\frac{\beta_2N_2^*-d_2-\gamma_2}{\beta_2+(d_2+\gamma_2)\alpha_2}| < \epsilon$;

if $I_2(T) > \frac{\beta_2N_2^*-d_2-\gamma_2}{\beta_2+(d_2+\gamma_2)\alpha_2}$, there are two cases:

ⅰ) $I_2(t)$ is decreasing for $t > T$. In this case, there exists $T^*_2 > T$, such that $|I_2(T^*_2)-\frac{\beta_2N_2^*-d_2-\gamma_2}{\beta_2+(d_2+\gamma_2)\alpha_2}| < \epsilon$;

ⅱ) There exists $T_1 > T$, such that $\frac{dI_2(T_1)}{dt} > 0$. That is

Since $I_1(t) < \epsilon_1$ for $t > T$, we have

Since $I_2(T) > \frac{\beta_2N_2^*-d_2-\gamma_2}{\beta_2+(d_2+\gamma_2)\alpha_2}$, we have

So

If only $\epsilon_1 < \frac{(d_2+\gamma_2)\alpha_2+\beta_2}{\left(\alpha_2+\frac{\beta_2+(d_2+\gamma_2)\alpha_2}{\beta_2N_2^*-d_2-\gamma_2}\right)m}\epsilon$, then $|I_2(T^*_2)-\frac{\beta_2N_2^*-d_2-\gamma_2}{\beta_2+(d_2+\gamma_2)\alpha_2}| < \epsilon.$ It is completed.

3.4.   The stability of the epidemic equilibrium

In this subsection, we will discuss the local stability of the epidemic equilibrium firstly and then discuss the global asymptotical stability.

Theorem 3.6. For the system (2.5), if $R_0^1 > 1$, the epidemic equilibrium $E_{*}$ is locally asymptotically stable.

Proof. The Jacobian matrix at the epidemic equilibrium $E_{*}$ of system (2.5) is

The two eigenvalues of the Jacobian matrix are

and

It is easy to see that if $R_0^1 > 1$, $\lambda_1 < 0$. The next, we need only prove the second eigenvalue $\lambda_2 < 0$ if $R_0^1 > 1$. Since $(1+\alpha_2I_2^*)^2 > 0$, we need only prove $\beta_2N_2^*-2\beta_2I_2^{*}-2\beta_2\alpha_2I_2^{*^2}-(d_2+\gamma_2)(1+\alpha_2I_2^{*})^2 < 0$. Let

Since $I_2^*\geq\frac{(\beta_2 N_2^*-d_2-\gamma_2+m\alpha_1I_1^*)}{(d_2+\gamma_2)\alpha_2+\beta_2}$, so

This completes the proof.

Theorem 3.7. For the systeml (2.5), if $R_0^1 > 1$, the epidemic equilibrium $E_{*}$ is globally asymptotically stable.

Proof. Since $E_{*}$ is stable when $R_0^1 > 1$, we need only prove $E_{*}$ is globally attractive.

Consider the equation

Let

Then $f_1^{\prime}(I_1) = \beta_1\frac{-1-N_1^*\alpha_1}{(1+\alpha_1I_1)^2} < 0$. So $f_1(I_1)$ is a monotonic decreasing function for all $I_1 > 0$. Furthermore, $f_1(0) > 0$, $f_1(N_1^*) < 0$ and $f_1(I_1^*) = 0$ when $R_0^1 > 1$. That means if $I_1\in(0, I_1^*),$ $f_1(I_1) > 0$ and $\dfrac{dI_1(t)}{dt} > 0$; if $I_1\in(I_1^*, N_1^*),$ $f_1(I_1) < 0$ and $\dfrac{dI_1(t)}{dt} < 0$. Hence $\lim\limits_{t\rightarrow \infty}I_1(t) = I_1^*$. By Eq (3.3), $\lim\limits_{t\rightarrow \infty}I_2(t) = I_2^*$. Thus $E_{*}$ is globally asymptotically stable.

The results about the existence and stability of equilibria are summarized in Table 1.

Remark 3.3. From Theorems 3.5 and 3.7, we know that for the system (2.5), if $R_0 > 1$, the infectious disease is uniformly persistent. However, the infectious disease is not always uniformly persistent in every patch. If $R_0 > 1$, but $R_0^1 < 1$, the disease is uniformly persistent in the sink patch, but is extinct in the source patch. If $R_0^1 > 1$, the disease is always uniformly persistent in every patch. This is a different conclusion resulted by the reducible migration matrix.

Remark 3.4. From Theorem 3.5, in the case that $R_0^1 < 1$ and $R_0^2 > 1$, the infection does not persist in the source patch but is able to persist in the sink patch. So, in the early stage of the spread of infectious disease, the sink patch should assess the reproduction number $R_0^2$ reasonably and take control measures timely to prevent the epidemic.

4.   Simulations

In this section, we carry on numerical simulations to verify the theoretical conclusions, reveal the influence of the migration rate form patch 1 to patch 2 on the basic reproduction number, the transmission scale and transmission speed, and discuss the influence of the parameters $\alpha_1$ and $\alpha_2$ that measure the inhibitory effect on the basic reproduction number, the transmission scale and transmission speed.

To numerically illustrate the theoretical results, we need to choose some parameter values (see Table 2).

We verify the theoretical conclusions firstly. Denote the density of the infective individuals in patch 1 by $i_1(t) = \frac{I_1(t)}{N_1^*}$. Denote the density of the infective individuals in patch 2 by $i_2(t) = \frac{I_2(t)}{N_2^*}$. Figure 1 shows the evolution of the density of infective individuals in the two patches when $R_0^1 = 0.8889$ and $R_0^2 = 0.7500$. As predicted by the analytic calculation, the infectious disease in the two patches will disappear eventually. Figure 2 shows the evolution of the density of infective individuals in the two patches when $R_0^1 = 0.8889$ and $R_0^2 = 1.8750$. We can see the infectious disease will be endemic in patch 2 and the infectious disease in patch 1 will disappear eventually. Figure 3 shows the evolution of the density of infective individuals in the two patches when $R_0^1 = 1.4815$ and $R_0^2 = 2.2917$. We can see the infectious disease will be endemic in the two patches. And we can see the infectious disease will be endemic in the two patches if $R_0^1 > 1$ from the subfigures (c) and (d) of Figure 4.

Second, we reveal the influence of the migration rate $m$ on the transmission in Figure 5. With the increasing of $m$, the density of infective individuals in patch 1 $i_1$ is decreasing, however the density of infective individuals in patch 2 $i_2$ is increasing.

Third, we reveal the parameters $\alpha_1$ and $\alpha_2$ on the transmission scale and transmission speed. We can see that when $\alpha_1$ is increasing, the density of infective individuals in patch 1 $i_1$ is decreasing from Figure 6 and When $\alpha_2$ is increasing, the density of infective individuals in patch 2 $i_1$ is decreasing from Figure 7.

5.   Discussion and conclusions

Many scholars have studied infectious disease transmission with population migration ^{[1,2,3,4,5,6,7]}, assuming that the migration matrix is irreducible, and found that the propagation dynamics of infectious diseases is determined by the basic reproduction number of the system. When the basic reproduction number is less than 1, the infectious disease eventually becomes extinct; when the basic reproduction number is larger than 1, the infectious disease is epidemic eventually. Since the migration matrix is irreducible, all patches are a connected whole. In all patches, infectious diseases are either extinct or epidemic. That is there is not the phenomenon that infectious diseases are extinct in some patches but epidemic in the others.

Because the studies about the spread of infectious diseases with reducible migration matrix are rare, in this paper, we proposed a two-patch SIS model with saturating contact rate and one-directing population dispersal, discussed the global asymptotic stability of the disease-free equilibrium, the boundary equilibrium and the endemic equilibrium respectively, and revealed the influence of saturating contact rate and migration rate on basic reproduction number and the transmission scale. We have the following main conclusions:

1) If $R_0^1 > 1$ then the system tends to a global endemic equilibrium in which infected individuals are present in both patches provided initially there were infected individuals in the source patch; If $R_0^1 < 1$ and $R_0^2 > 1$ then the system converges to an equilibrium with infected individuals only in the sink patch; If $R_0^1 < 1$ and $R_0^2 < 1$ then the system converges to the disease-free equilibrium.

2) When migration rate is increasing, the density of infective individuals in the source patch is decreasing; but the density of infective individuals in the sink patch is increasing;

3) With the increasing of the parameter $\alpha_i$ $(i = 1, 2)$ in saturating contact rate, the density of infective individuals in patch $i$ $(i = 1, 2)$ is decreasing.

The similar conclusions can be obtained for the two patch SI model

We can generalize the current model in many aspects to increase realism. For instance, the infection rate can be given by $\beta(I)SI$. We can give the properties on function $\beta(I)$ such that $\beta(I)$ is decreasing and tends to 0 when I tends to infinity. The mortality rates of the susceptible and infected individuals are the same in the current model. In fact, the disease-induced death rate can not be neglected sometimes. So the disease-induced death rate can be considered. It is also significant to consider heterogeneous number of contacts for each individual on complex network. There are many paper on this topic ^[22,23]. One can investigate the multi-patch epidemic model with reducible migration matrix.

Acknowledgments

This work is supported by the National Natural Sciences Foundation of China (Nos.12001501, 12071445, 11571324, 61603351), Shanxi Province Science Foundation for Youths (201901D211216), the Fund for Shanxi '1331KIRT'.

Conflict of interest

All authors declare no conflicts of interest in this paper.