The backward bifurcation of an age-structured cholera transmission model with saturation incidence

1.   Introduction

Cholera patients have a large number of curved bacteria in their intestines. These bacteria are the pathogen of cholera, namely Vibrio cholerae. In other words, cholera is a clinical epidemiological syndrome caused by Vibrio cholerae ^[1]. There are two ways of transmission of cholera, environment-to-human and human-to-human ^[2,3]. For cholera control, the most important prevention strategy is based on traditional basic cholera prevention methods, such as improving health systems, safer water treatment and improving food and personal hygiene ^[4,5]. In addition, the World Health Organization recommends that an oral anti-cholera vaccine can be used in areas at risk of epidemic and outbreak ^[6,7,8].

Many mathematical models about cholera have been put forward to study the mechanism of cholera transmission and obtain better control measures ^[3,9,10,11]. In 2006, Hartley et al. brought Vibrio cholerae into a cholera model and got a better fitting result ^[3]. Based on the above research, in ^[9], the direct and indirect transmissions of cholera models were considered, and the importance of multiple transmission routes was obtained. In addition, some mathematical models considering age structures have been proposed ^[10,11]. In 2013, Brauer et al. proposed an age-structured cholera model, which considered the infection age of infected individuals and the biological age of pathogens ^[10]. In ^[12], Posny et al. put forward a cholera model with control measures, including basic hygiene, treatment and vaccine control strategies. Through analysis, it can be concluded that vaccines have a significant impact on disease control. However, the vaccine effect will decline with the increase of vaccination age. Therefore, the vaccination age of vaccinated individuals was incorporated into cholera transmission models ^[13,14]. In addition, Wang et al. proposed an age structured cholera model with multiple transmission routes, in which the infection age of the infected person and the biological age of the pathogen in the environment are considered ^[15].

In ^[13,14], the authors assume that the incidence rate is bilinear, that is, the infection term is a linear increasing function of the number of infected individuals and the concentration of V. cholerae. However, for human-to-human transmission, as the number of infected individuals increases, the number of contacts of a susceptible per unit time cannot always increase linearly with the number of infected individuals. Similarly, for environment-to-human transmission, the concentration of Vibrio cholerae in the environment is gradually saturated. Therefore, in ^[16], Capasso and Serio introduced a saturation incidence rate $\frac{\beta I}{1+\alpha I}$ to measure the crowding effect. In addition, in ^[17], Codeco considered the saturation incidence for environment-to-human transmission and studied the impact of saturation incidence on cholera transmission.

Based on the above work, in this paper, we consider a cholera model with saturation incidence, the vaccination age of vaccinated individuals, the infection age of infected individuals and the biological age of pathogens. In addition, we assume that the recovered individuals can be infected again. Here, the reason why we consider the biological age of the pathogen in our model is that the Vibrio cholerae shed from the human gastrointestinal tract has been proved by experiments to be highly infectious, and the hyper-infectious state was transient ^[3]. With the passage of time, the infectivity will gradually decrease, from high infectivity to low infectivity.

The structure of this paper is organized as follows. In Section 2, we propose an age-structured cholera model with saturation incidence, vaccination age of vaccinated individuals, infection age of infected individuals and biological age of pathogens. In addition, the basic properties on the positivity and boundness of solutions are discussed. In Section 3, we get the asymptotic smoothness of the semi-flow of the system. In Section 4, we calculate the basic reproduction number and study the local stability of the disease-free equilibrium when $\mathcal{R}_{0} < 1$. In Section 5, we obtain the existence of backward bifurcation of the model. In Section 6, we carry out numerical simulations to illustrate our analytical results and study the impact of various control measures on the number of infected individuals. The paper ends with a brief discussion in the last section.

2.   Model formulation

Let $S(t)$ denote the density of susceptible individuals at time $t$, $V(a, t)$ denote the density of the vaccinated individuals with vaccination age $a$ at time t, $i(b, t)$ denote the density of the infected individuals with infection age $b$ at time $t$, $R(t)$ denote the density of recovered individuals at time $t$ and $B(c, t)$ denote the concentration of V. cholera with biological age $c$ at time $t$. We assume that susceptible individuals are recruited at a constant rate $A$, and susceptible individuals can acquire cholera infection by environment-to-human transmission at a rate $\int_{0}^{\infty}\frac{\beta_{1}(c)B(c, t)}{K+B(c, t)}dc$ and human-to-human transmission at a rate $\int_{0}^{\infty}\frac{\beta(b)i(b, t)}{1+\alpha i(b, t)}db$. The vaccination rate of susceptible individuals is $\phi$, and all human individuals have a natural death rate $\mu$. For vaccinated individuals, there is the reduction of vaccine efficacy $\gamma_{1}(a)$ depending on the vaccination age $a$. Infected individuals have an age-dependent removal rate $\theta(b) = \mu+\gamma_{2}(b)+\rho(b)$, where $\gamma_{2}(b)$ is the recovery rate and $\rho(b)$ is the disease-induced death rate. The acquired immunity of recovered individuals is assumed to wane at rate $\omega$. For the environmental bacteria V. cholera, there are age-dependent shedding rate $\xi(b)$ and removal rate $\delta(c)$. In addition, it should be noted that $V(t) = \int_{0}^{\infty}V(a, t)da, I(t) = \int_{0}^{\infty}i(b, t)db, B(t) = \int_{0}^{\infty}B(c, t)dc.$

Based on these assumptions, in this paper, we consider the following differential equations:

with the boundary conditions

and initial conditions

Assumption 1.1 The parameters satisfy the following conditions.

$\bf 1)$ $\gamma_{1}(a)$, $\beta(b)$, $\beta_{1}(c)$, $\theta(b)$, $\gamma_{2}(b)$, $\xi(b)$, $\delta(c)\in L_{+}^{1}(0, \infty)$ with positive essential upper bounds $\bar{\gamma_{1}}$, $\bar{\beta}$, $\bar{\beta_{1}}$, $\bar{\theta}$, $\bar{\gamma_{2}}$, $\bar{\xi}$, $\bar{\delta}$, respectively.

$\bf 2)$ $\beta(b)$ and $\beta_{1}(c)$ are bounded, uniformly continuous and with compact support.

Denote the function space $X_{0} = R\times L^{1}(0, \infty) \times L^{1}(0, \infty)\times R\times L^{1}(0, \infty)$, equipped with the norm

In addition, we define the positive cone $X_{0}^{+} = R_{+}\times L_{+}^{1}(0, \infty) \times L_{+}^{1}(0, \infty)\times R_{+}\times L_{+}^{1}(0, \infty)$.

Lemma 2.1. For every $\phi_{0} = (S_{0}, V_{0}(a), i_{0}(b), R_{0}, B_{0}(c))\in X_{0}^{+}$, model $(2.1)$ has a unique nonnegative solution.

Proof. Let $\Omega$ be any bounded set in $X_{0}^{+}$. Then, there exists a positive number $M$ such that for any $\phi_{0}\in\Omega, ||\phi_{0}||_{X_{0}^{+}}\leq M.$ First, we define a nonlinear operator $G(S(t), R(t), \textbf{u}(t))$ by

where $G_{i}$ is defined on $R^{2}\times(L^{1}(0, \infty))^{3}, i = 1, 2.$

In addition, we define

and

where $B, F$ are defined on $R^{2}\times(L^{1}(0, \infty))^{3}$. Based on the above definitions, we can translate model (2.1) into the following abstract form:

Therefore, for any $\phi_{0}, \bar{\phi}_{0}\in\Omega$,

where $L_{G} = \max\{\mu+\phi+\bar{\beta}M+\frac{\bar{\beta}_{1}M}{K}, \omega, \bar{\gamma}_{1}\}.$ Similarly, $B[S(t), \textbf{u}(t)]$ and $G_{2}[R(t), \textbf{u}(t)]$ satisfy

where $L_{B}$ and $L_{G_{2}}$ are Lipschitz coefficients.

Furthermore, we have

and $B[\phi(t)]\geq0$ for all $\phi_{0}\in X_{0}^{+}.$ Hence, Conditions 1–5 in Theorem 2.1 ^[18] are satisfied and model (2.1) has a unique nonnegative solution.

Next, for the sake of convenience, define

Solving model (2.1) along the characteristic lines $t-a =$ constant, $t-b =$ constant and $t-c =$ constant, yields

and

respectively.

Theorem 2.1. All solutions of model $(2.1)$ with boundary conditions $(2.2)$ and initial conditions $(2.3)$ are ultimately bounded.

Proof. Let $N(t) = S(t)+\int_{0}^{\infty}V(a, t)da+\int_{0}^{\infty}i(b, t)db+R(t)$. We obtain that

It follows from (2.5) that

Since $\pi_{1}(0) = 1$ and $d\pi_{1}(a)/da = -(\mu+\gamma_{1}(a))\pi_{1}(a)$, we have

Similarly, we have

Therefore, we can get

which yields

Hence, for $\varepsilon > 0$ sufficiently small, there is a $T_{1} > 0$ such that if $t > T_{1}$,

Similarly, we have

which yields

Thus this completes the proof.

Thus, we can define the solution semi-flow $\Phi(t, x_{0}) = \Phi_{t}(x_{0}) = (S(t), V(\cdot, t), i(\cdot, t), R(t), B(\cdot, t)), t\geq0, x_{0} = (S_{0}, V_{0}(a), i_{0}(b), R_{0}, B_{0}(c))$. The functions $V_{0}(a), i_{0}(b), B_{0}(c)$ are in $L_{+}^{1}(0, \infty).$

3.   Asymptotic smoothness

In this section, we investigate the asymptotic smoothness of the semi-flow $\Phi(t, x_{0})$.

Lemma 3.1. $S(t), f_{1}(t)$ and $f_{2}(t)$ are Lipschitz continuous on $R^{+}$ with Lipschitz coefficients $L_{S}, L_{1}$ and $L_{2}$, respectively, that is, the following inequalities hold:

where $L_{1} = [(\bar{\beta}+\bar{\beta}_{1}/K)M+\bar{\beta}\bar{\theta}+L_{\beta}]M, L_{2} = \frac{1}{K}(\bar{\beta}_{1}\bar{\xi}+\bar{\beta}_{1}\bar{\delta}+L_{\beta_{1}})M,$ in which $L_{\beta}$ and $L_{\beta_{1}}$ are the Lipschitz coefficients of $\beta$ and $\beta_{1},$ respectively.

Theorem 3.1. The semi-flow ${\Phi(t, x_{0})}$ generated by model $(2.1)$ is asymptotically smooth.

Proof. First, the semi-flow ${\Phi(t, x_{0})}$ is expressed as the sum of the following two operators $\phi(t, x_{0}), \varphi(t, x_{0}):R^{+}\times X_{0}\rightarrow X_{0},$ where

with

For $h > 0$, let $u(t, h) = he^{-\mu_{0}t}$. It is obvious that $\lim\limits_{t\rightarrow +\infty}u(t, h) = 0$. We obtain that

that is,

If $||x_{0}||_{X_{0}} < h$, then $||\phi_{t}(x_{0})||_{X_{0}}\leq he^{-\mu_{0}t}\triangleq u(t, h)$. Thus, the condition $(i)$ of Lemma 3.2.3 in ^[19] is satisfied. On the other hand,

Note that

Hence, we can obtain that

In addition,

According to Lemma 3.1,

It follows from (3.1), (3.2) and (3.3)that

Hence, we can conclude that $V_{2}(a, t)$ satisfies the conditions of Lemma 3.2.3 in ^[19]. Similarly, we also conclude that $i_{2}(b, t)$ and $B_{2}(c, t)$ satisfy the conditions of Lemma 3.2.3 in ^[19]. According to Theorem 27 in ^[20], the $\varphi(t, x_{0})$ is completely continuous. Therefore, the semi-flow ${\Phi(t, x_{0})}$ is asymptotically smooth.

4.   Equilibria analysis

It is easy to show that model (2.1) has a disease-free equilibrium $E^{0}(S^{0}, V^{0}(a), 0, 0, 0),$ where

where $D(0) = \int_{0}^{\infty}\gamma_{1}(a)\pi_{1}(a)da$, and we have

Linearizing model (2.1) at disease-free equilibrium $E^{0}$, we obtain

with the boundary conditions

Further, we have

where

Equation (4.3) generates a classical renewal equation. The basic reproduction number is defined as

where $m = \int_{0}^{\infty}\xi(b)\pi_{2}(b)db.$ The first term $S^{0}\int_{0}^{\infty}\beta(b)\pi_{2}(b)db$ refers to the number of secondary infections caused by infected individuals during the infection period, and the second term $S^{0}\frac{m}{K}\int_{0}^{\infty}\beta_{1}(c)\pi_{3}(c)dc$ refers to the number of secondary infections caused by Vibrio cholera exposed to the environment by susceptible individuals.

Theorem 4.1. The disease-free equilibrium $E_{0}$ is locally asymptotically stable when $\mathcal{R}_{0} < 1$, and unstable when $\mathcal{R}_{0} > 1$.

Proof. Letting $x_{1}(t) = S(t)-S^{0}, x_{2}(a, t) = V(a, t)-V^{0}(a), x_{3}(b, t) = i(b, t), x_{4}(t) = R(t), x_{5}(c, t) = B(c, t)$ and linearizing model (2.1) at $E_{0}$, we obtain

Substituting $x_{1}(t) = \widetilde{x}_{1}e^{\lambda t}, x_{2}(a, t) = \widetilde{x}_{2}(a)e^{\lambda t}, x_{3}(b, t) = \widetilde{x}_{3}(b)e^{\lambda t}, x_{4}(t) = \widetilde{x}_{4}e^{\lambda t}, x_{5}(c, t) = \widetilde{x}_{5}(c)e^{\lambda t}$ into (4.5), we have

From $(4.6a)$ and $(4.6g)$, we have $\widetilde{x}_{3}(0) = -(\lambda+\mu+\phi)\widetilde{x}_{1}+\omega\widetilde{x}_{4}+\int_{0}^{\infty}\gamma_{1}(a)\widetilde{x}_{2}(a)da$. Then, integrating (4.6c) from $0$ to $b$ yields

According to (4.6e), integrating (4.6e) from $0$ to $c$ yields

Substituting (4.7) and (4.8) into $(4.6g)$, we obtain that

Dividing $\widetilde{x}_{3}(0)$ from both sides of Eq (4.9), we have

Next, we claim that all roots of (4.10) have negative real parts if $\mathcal{R}_{0} < 1.$ Otherwise, there exists a root $\lambda_{1} = x_{1}+iy_{1}$ with $x_{1}\geq0$. In this case, substituting $\lambda_{1}$ into (4.10), we obtain

It follows that

which contradicts with (4.11). Thus, if $\mathcal{R}_{0} < 1$, $E_{0}$ is locally asymptotically stable.

If $\mathcal{R}_{0} > 1$, let

It is easy to obtain for the expression $H(\lambda)$ that

From the above, we know that if $\mathcal{R}_{0} > 1$, $H(\lambda)$ has at least a root with positive parts. Therefore, $E_{0}$ is unstable.

5.   The existence of backward bifurcation

Theorem 5.1. Model $(2.1)$ has at least one endemic equilibrium $E^{*}(S^{*}, V^{*}(a), i^{*}(b), R^{*}, B^{*}(c))$ when $\mathcal{R}_{0} > 1.$

Proof. Any endemic equilibrium $E^{*}(S^{*}, V^{*}(a), i^{*}(b), R^{*}, B^{*}(c))$ of model (2.1) satisfies the following equations:

From (5.1), we obtain

Substituting (5.2) into $i^{*}(0) = S^{*}\left(\int_{0}^{\infty}\frac{\beta(b)i^{*}(b)}{1+\alpha i^{*}(b)}db+\int_{0}^{\infty}\frac{\beta_{1}(c)B^{*}(c)}{K+B^{*}(c)}dc\right)$, we obtain

Dividing $i^{*}(0)$ from both sides of this equation, we have

where

in which

$N(i^{*}(0)) = \left(\int_0^{\infty}\frac{\beta(b)\pi_{2}(b)}{1+\alpha i^{*}(0)\pi_{2}(b)}db+\int_0^{\infty}\frac{\beta_1(c)m\pi_{3}(c)}{K+i^{*}(0)m\pi_{3}(c)}dc\right)$ and $m = \int_{0}^{\infty}\xi(b)\pi_{2}(b)db.$ Let $g(i^{*}(0)) = S^{*}N(i^{*}(0))-1,$ where $g(0) = S^{0}N(0))-1 = \mathcal{R}_{0}-1$, and $g(\infty) < 0.$ It is easy to show that if $\mathcal{R}_{0} > 1$, we have $g(S^{0}) > 0$. Thus, (5.4) has at least one positive root. Model (2.1) has at least one endemic equilibrium $E^{*}(S^{*}, V^{*}(a), i^{*}(b), R^{*}, B^{*}(c))$ when $\mathcal{R}_{0} > 1.$

Letting $y_{1}(t) = S(t)-S^{*}, y_{2}(a, t)$$= V(a, t)-V^{*}(a),$$y_{3}(b, t) = i(b, t)-i^{*}(b),$$y_{4}(t) = R(t)-R^{*},$$y_{5}(c, t) = B(c, t)-B^{*}(c)$ and linearizing model (2.1) at $E^{*}$, we obtain

Substituting $y_{1}(t) = \widetilde{y}_{1}e^{\lambda t}, y_{2}(a, t) = \widetilde{y}_{2}(a)e^{\lambda t}, y_{3}(b, t) = \widetilde{y}_{3}(b)e^{\lambda t}, y_{4}(t) = \widetilde{y}_{4}e^{\lambda t}, y_{5}(c, t) = \widetilde{y}_{5}(c)e^{\lambda t}$ into (5.5), we have

From $(5.6a)$ and $(5.6g)$, we have $\widetilde{y}_{3}(0) = -(\lambda+\mu+\phi)\widetilde{y}_{1}+\omega\widetilde{y}_{4}+\phi \tilde{y}_{1}\int_{0}^{\infty}\gamma_{1}(a)e^{-a\lambda}\pi_{1}(a) \text{d}a$. Then, integrating (5.6c) from $0$ to $b$ yields

According to (5.6e), integrating (5.6e) from $0$ to $c$ yields

Substituting (5.7) and (5.8) into $(5.6g)$, we obtain that

If all roots of the characteristic Eq (5.9) have negative real parts, the endemic equilibrium $E^{*}$ is locally asymptotically stable, otherwise, $E^{*}$ is unstable.

Now, we study the existence of backward bifurcation of model (2.1). We derive the necessary and sufficient conditions of backward bifurcation in model (2.1) by using the method of the sign of the derivative. In order to apply the sign of derivative method, first, we need to select a bifurcation parameter and key infectious quantity. Here, we choose $\bar{\beta}$ as a bifurcation parameter and $i(0)$ as a key infectious quantity. In addition, we rewrite $\beta(b) = \bar{\beta}\beta_{0}(b)$ where $\sup\beta_{0}(b) = 1$. From the equation of $i^{*}(0)$, we have

where

Differentiating with respect $\bar{\beta}$ and calculating at $E_{0}$, we have

where

Hence, $i'(0) < 0$ if and only if

Therefore, we obtain that the necessary and sufficient conditions for the backward bifurcation of model (2.1). The model (2.1) exhibits backward bifurcation if there is a set of parameters such that (5.10) is satisfied. Clearly, since the integral is positive we can choose $A, \alpha$ small enough and $q$ big enough so that the inequality holds.

6.   Numerical results

In this section, we perform some numerical simulations to illustrate our analytical results. We also compare the effects of various control measures, including basic control measures and vaccination on the number of infected individuals. All parameter values are defined in (6.1) and Table 1 except the vaccination rate $\phi$.

Based on the reduction of vaccine efficacy $\gamma_{1}(a)$, the recovery rate $\gamma_{2}(b)$, the transmission coefficient $\beta(b)$, the removal rate of infected individuals $\theta(b)$ and the removal rate of V. cholera $\delta(c)$ are increasing functions; the transmission coefficient $\beta_{1}(c)$ and the shedding rate of V. cholera $\xi(b)$ are decreasing functions. Here, we assume $\gamma_{1}(a), \gamma_{2}(b), \beta(b), \beta_{1}(c), \theta(b), \xi(b), \delta(c)$ satisfy the following expressions:

Taking $\gamma_{1}(a)$ as an example, 0.0530 is the initial value of reduction of vaccine efficacy $\gamma_{1}(a)$, 0.0866 is the final value of reduction of vaccine efficacy $\gamma_{1}(a)$, and 0.8564 is the exponential rate.

Through theoretical analysis, we obtain that the disease-free equilibrium is locally asymptotically stable if $\mathcal{R}_{0} < 1$. Here, we choose $\phi = 0.1$. By calculation, we get the basic reproduction number $\mathcal{R}_{0} = 0.8425 < 1$. As shown in Figure 1, when $\mathcal{R}_{0} < 1$, the number of infected individuals $I(t) = \int_{0}^{\infty}i(b, t)db$ and the concentration of Vibrio cholerae in the environment $B(t) = \int_{0}^{\infty}B(c, t)dc$ converge to the level at the disease-free equilibrium $E_{0}$. In addition, if we choose $\phi = 0.008$, by calculation, we get $\mathcal{R}_{0} = 2.9458 > 1$. As shown in Figure 2, when $\mathcal{R}_{0} > 1$, the number of infected individuals $I(t) = \int_{0}^{\infty}i(b, t)db$ and the concentration of Vibrio cholerae in the environment $B(t) = \int_{0}^{\infty}B(c, t)dc$ converge to the respective levels at the endemic equilibrium.

Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficient (PRCC) are usually adopted to identify the uncertainties of model parameters, and we can use PRCC to evaluate the influence of parameters in the model on the basic reproduction number $\mathcal{R}_{0}$. Here, for $\gamma_{1}(a), \beta(b),$ $\beta_{1}(c),$ $\theta(b)$, $\xi(b)$, $\delta(c)$, we take their values as constants (the average of the initial and final values). From Figure 3, we see that the vaccination rate $\phi$ has a great impact on the basic reproduction number. This indicates that the measures related to the vaccine have a great impact on the basic reproduction number $\mathcal{R}_{0}$.

Different control strategies have different effects on the number of infected individuals. First, we consider the impact of the intensity of basic control measures on the number of infected individuals. Here, we choose $\phi = 0.005$. Among the parameters of the model, the basic anti-cholera control measures are determined by the parameters $\beta(b), \beta_{1}(c)$. Reductions in $\beta(b), \beta_{1}(c)$ imply reductions in the rate of human exposure to cholera pollution sources and the rate of human contact with infected individuals, respectively. This can be reduced by strengthening some basic control measures, that is, the smaller the values of $\beta(b), \beta_{1}(c)$, the higher the degree of basic control measures, and the higher the value of $\beta(b), \beta_{1}(c)$, the lower the degree of basic control measures. Hence, we select different parameter values to analyze the impact of basic control measures on the number of infected individuals $I(t) = \int_{0}^{\infty}i(b, t)db$ (see Figure 4).

Next, we consider the effect of vaccination on model dynamics. We analyze the dynamic behavior of the model from the aspect of vaccination rate. We select different vaccination rates $\phi$ and compare their effects on the number of infected individuals $I(t) = \int_{0}^{\infty}i(b, t)db$. From Figure 5, we can see that the increase of the vaccination rate can reduce not only the total number of infected individuals but also the duration of the epidemic.

7.   Discussion

In this paper, we have presented a cholera model with the vaccination age of vaccinated individuals, the infection age of infected individuals and the biological age of pathogen, and we investigated the dynamics of cholera transmission. By calculation, we obtain the basic reproduction number. When $\mathcal{R}_{0} < 1$, the local stability of the disease-free equilibrium is obtained. In addition, we also obtain the existence of backward bifurcation of system (2.1). In the numerical simulations, we verify the results of our theoretical analysis (see Figure 1) and analyze the impact of various control measures on the number of infected individuals. Firstly, we consider the impact of basic control measures on infected individuals (see Figure 4). From Figure 4, we can see that the basic control measures can reduce the number of infected individuals to a certain level. However, it is difficult to achieve the effect of clearing the number of infected individuals. According to the sensitivity analysis, we can observe that the vaccination rate $\phi$ has a great influence on the basic reproduction number $\mathcal{R}_{0}$ (see Figure 3). That is to say, the vaccination rate has an important influence on the cholera transmission control. Further, Figure 5 analyzes the impact of the vaccination rate $\phi$ on the number of infected individuals. From Figure 5, we can see that increasing the vaccination rate can achieve the good effect of clearing the number of infected individuals.

It should be pointed out that the model we are considering here is the classic epidemic model. However, the classical epidemic model assumes that the population is homogeneously mixed, and it is more realistic to consider the heterogeneous number of contacts of each individual, that is, the network epidemic model ^[21,22]. Therefore, extending the model to complex networks can be more realistic.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (No. 61873154; 11331009), Shanxi Key Laboratory (No. 201705D111006), Shanxi Scientific and Technology Innovation Team (No. 201805D131012-1).

Conflict of interest

All authors declare no conflicts of interest in this paper.