Dynamical analysis and optimal control of an multi-age-structured vector-borne disease model with multiple transmission pathways

1.   Introduction

According to a report by the World Health Organization (WHO), vector-borne diseases account for approximately 17% of all illnesses caused by infectious diseases. Common vector-borne diseases are Malaria, Dengue fever, Zika, Cholera, Yellow fever, West Nile fever, and so on. Vector-borne diseases result in more than 1 billion illnesses and over 1 million fatalities annually ^[1]. Considering that vector-borne diseases have caused great damage to human health and the social economy, it is urgent to study the laws of their transmission and how to prevent and control them.

In the past few decades, many scholars have conducted in-depth and detailed studies on the spread and control of various vector-borne diseases by establishing mathematical models and a series of remarkable research results have been achieved ^[2,3,4,5,6]. Major studies include the expression of the basic reproduction number, the existence and stability of various equilibria, the persistence and extinction of diseases, and optimal control problems, etc. In particular, Brown et al.^[7] presented a mathematical model of Cholera transmission, established the local and global stability of the equilibria characterizing the threshold dynamics of Cholera, and discussed the influence of pathogens in the dynamics of Cholera spread. In^[8], Kuddus et al. proposed a model for the dynamics of Malaria transmission between humans and mosquitoes, analyzed the existence and global stability of equilibria and found that the rate of contact between humans and mosquitoes had a significant effect on Malaria spread through parameter estimation. Wang et al. ^[9] studied the threshold dynamics of a Dengue epidemic model, and the results showed that controlling the number and activity of vectors can significantly affect the speed and extent of disease spread. Further, it is noted that vector-borne diseases can be transmitted not only through vectors but also between individuals and individuals. For example, Chen et al. ^[10] found that Zika can be transmitted not only between hosts and vectors to each other but also within host populations. More specifically, an infected person can spread it to his (her) partner through sexual intercourse.

It is well known that the physiological age of the host population is a key factor in the transmission and control of infectious diseases, as the risk of infection varies with age, and interactions between different age groups are heterogeneous. For example, SARS-CoV-2, the virus behind COVID-19, is particularly dangerous for individuals over 60, especially those over 80. Centers for Disease Control reports show that 31–59% of 75–84-year-olds diagnosed require hospitalization, compared to 14–21% of those aged 20–44, highlighting the influence of age in disease severity and transmission dynamics. Noting this feature, researchers have proposed various age-structured models to study the transmission dynamics of infectious diseases ^[11,12,13]. In ^[14], Cai et al. established an age-structured Cholera model and discussed the threshold dynamic behavior of the model and the effect of actual age structure on the spread trend of the disease. Huang et al. ^[15] presented a stability study of an age-structured epidemiological model and showed that actual age affects an individual's risk of infection and more accurately describes the dynamics of disease transmission in different age groups. Yu et al. ^[16] proposed an age-structured COVID-19 model, discussed the significant differences in the influence of different age groups on the spread of disease due to differences in their activity ranges, and revealed the importance of age structure in the disease transmission process.

It is worth noting that in the process of transmission of vector-borne diseases, the spread capacity of the infected vector is closely related to its age of infection. Particularly for vectors such as mosquitoes with short life cycles, which have different transmission potentials at different stages of post-infection. Therefore, in order to analyze the dynamics of disease spread more accurately, it is essential to consider that mosquitoes have an age of infection. For example, Liang et al. ^[17] presented a model of vector-borne disease with multiple class age structures, gave an exact expression for the basic reproduction number, and discussed the effect of the age of infection of the vectors on the basic reproduction number and the spread and control of the disease. In ^[18], Richard et al. proposed a model of human-vector Malaria transmission related to infection age, discussed the effect of vector infection age on Malaria transmission, and found through numerical simulations that the effects of different intervention policies differed in mosquito populations with different ages of infection.

The issue of optimal control is of great significance in the development of preventive measures for the spread of disease. Therefore, the study of the optimal control problem for age-structured epidemiological modeling has attracted the attention of many researchers ^{[19,20,21,22]}. For example, Lin et al. ^[23] proposed an age-structured Cholera model with vaccination as a control strategy; The results of the study suggest that vaccination strategy at the beginning of a Cholera outbreak can significantly reduce the number of infections and is the most cost-effective strategy. Wang et al. ^[24] proposed and studied the global dynamics of an age-structured Malaria model with vaccination; The existence of optimal control is analyzed and effective measures to control Malaria transmission are obtained. Khan et al. ^[25] developed an age-structured SEIR model, using vaccination and treatment as control measures; The existence of optimal control variables was demonstrated using a suitable objective function.

Based on the above discussion, a novel coupled ordinary differential-partial differential equations model is proposed to discuss the impact of physiological age, infection age, multiple transmission routes, and asymptomatic infected persons on the spread of vector-borne diseases. The structure of our paper is arranged as follows. In Section 2, the model is given, and the existence and uniqueness of a global positive solution are verified. In Sections 3 and 4, the basic reproduction number $\mathcal{R}_0$ is defined, and the global asymptotic stability of the disease-free steady state and the local asymptotic stability of the endemic steady state are discussed. Section 5 demonstrates the consistent persistence of the model. In Section 6, vaccination and insecticide spraying are applied to the model to demonstrate the existence of optimal control. The main results are explained through numerical simulations in Section 7, and a short conclusion is given in the last section.

2.   Model formulation

The host population of an area is classified into four mutually exclusive categories: susceptible, asymptomatically infected, symptomatically infected, and recovered individuals. And, their density functions with age $a$ at time $t$ are denoted as $S_{h}(t, a)$, $A_{h}(t, a)$, $I_{h}(t, a)$, and $R_{h}(t, a)$, respectively. Therefore, the total host population is given by $N_{h}(t, a) = S_{h}(t, a)+A_{h}(t, a)+I_{h}(t, a)+R_{h}(t, a)$. The vector population is divided into two subclasses: susceptible and infected vectors, where the quantity of susceptible vectors at time $t$ is denoted by $S_v(t)$ and the density function of infectious vectors with infection age $b$ at time $t$ is denoted by $I_v(t, b)$. Then, the total quantity of vector population is given by $N_{v}(t) = S_{v}(t) + \int_{0}^{\infty}I_{v}(t, b)\mathrm{d}b$. According to the transmission law of pathogens between host populations and vectors, susceptible individuals can be infected by the bite of infected vectors at a rate of $\lambda_{1}(t, a)$, and susceptible host can also be infected by contacting the infected hosts at a rate $\lambda_{2}(t, a)$, where the force of infection is defined as follows:

here, $z_1(a)$ and $z_2(a)$ represent the age-dependent contact rates of vector-to-host and host-to-host, respectively; $\beta_1(b)$ and $\beta_2(a)$ denote the infection rates of vector-to-host and host-to-host, respectively. Further, considering the host behavioral habits, it is assumed that $z_1(a) = \rho z_2(a)$, where $\rho$ is a positive constant.

Based on the above statements, a novel epidemic model with multiple transmission routes and multiple age-factor couplings is constructed

with the boundary and initial conditions

The biological explanations for the other parameters of model (2.1) are shown in Table 1.

Remark 1. Due to the prevalence of vector-borne diseases, numerous scholars have established and discussed vector-borne diseases ^[18,26]. In our model, considering the behavioral differences of hosts of different ages, it is proposed that host populations have physiological ages. Given the relatively short life cycle of vector populations, it is not necessary to consider their physiological age. However, there is a strong correlation between infectivity and age at infection, so it is important to consider the age of infection of the vector.

In addition, the following assumptions are reasonable based on the biological background of (2.1).

($\mathrm{H}_1$) $\Lambda_v$ and $\mu_v$ are positive constants, $q\in(0, 1)$, and functions $S_{h0}(a)$, $A_{h0}(a)$, $I_{h0}(a)$, $R_{h0}(a)$, and $I_{v0}(b)$ are non-negative, continuously integrable functions.

($\mathrm{H}_2$) Functions $b_h(a)$, $\beta_1(a)$, $\beta_2(a)$, $\beta_3(a)\in L^{1}_{+}(\mathbb{R}_+)$ and they are extended to zero outside the maximum age $a^{\dagger}$, where, $L_{+}^{1}(\mathbb{R}_+)$ is the space of Lebesgue integrable nonnegative functions on the interval $[0, +\infty)$.

($\mathrm{H}_3$) Functions $k(a)$, $\gamma_1(a)$, $\gamma_2(a)$, $\mu(a)\in L^{1}_{+}(\mathbb{R}_+)$ and $\int_{0}^{\infty}\varphi(a)\mathrm{d}a = +\infty$, where, $\varphi(a) = k(a)$, $\gamma_1(a)$ and $\gamma_2(a)$.

($\mathrm{H}_4$) There is a positive constant $\mu_{h0}$ such that $\mu_h(a)\geqslant \mu_{h0}$ for $a\in[0, a^{\dagger}]$.

Notice that the total quantity of vector population satisfies

Therefore, $N_{v}(t) = N_{v}(0)\mathrm{e}^{-\mu_vt}+\frac{\Lambda_v}{\mu_v}(1-\mathrm{e}^{-\mu_vt})$, that is, $\lim_{t\rightarrow \infty}N_{v}(t) = \frac{\Lambda_v}{\mu_v}$. According to the limiting theory of dynamical systems ^[27], the dynamics of model (2.1) is equivalent to the following model:

Next, we consider the existence and uniqueness of the global nonnegative solutions of model (2.3). To do so, define a Banach space $\mathbb{X}^n = L^1(\mathbb{R}_+)\times \cdots\times L^1(\mathbb{R}_+)$ with the norm $\|\varphi\| = \sum_{i = 1}^{n}\|\varphi_i\|$ for $\varphi = (\varphi_1, \cdots, \varphi_n)\in\mathbb{X}^n$, where $\|\varphi_i\| = \int_{0}^{\infty}|\varphi_i(a)|\mathrm{d}a$. Obviously, $\mathbb{X}^n_+ = L_+^1(\mathbb{R}_+)\times \cdots \times L_+^1(\mathbb{R}_+)$ is the positive cone of $\mathbb{X}^n$. Further, we denote $\mathbb{Y}_1 = \mathbb{X}_+^2$ and the linear operator $\mathcal{A}_{v, 1}:D(\mathcal{A}_{v, 1})\subset \mathbb{Y}_1\to \mathbb{Y}_1$ be defined by

where the domain $D(\mathcal{A}_{v, 1})$ of operator $\mathcal{A}_{v, 1}$ is $D(\mathcal{A}_{v, 1}) = \{0_{L^1(\mathbb{R}_+)}\}\times W^{1, 1}(\mathbb{R}_+)$, $W^{1, 1}(\mathbb{R}_+)$ represents the Sobolev space of all absolutely continuous functions on $\mathbb{R}_+$.

According to the assumptions ($\mathrm{H}_1$), for $\lambda\in\mathbb{C}$ ($\mathbb{C}$ is the domain of complex numbers) with $Re (\lambda)\geq-\mu_v$, then $\lambda\in r(\mathcal{A}_{v, 1})$ ($r(\mathcal{A})$ is the resolvent of an operator $\mathcal{A})$. Therefore, the explicit expression for the resolvent of the operator $\mathcal{A}_{v, 1}$ is derived

for $(\psi_1, \psi_2)^T\in \mathbb{X}^2_+$, where $\mathbb{I}$ represents the unit matrix. Further, one defines linear operators $\mathcal{A}_{h, i}:D(\mathcal{A}_{h, i})\subset \mathbb{Y}_{1}\to \mathbb{Y}_{1}$, $i = 1, \ 2, \ 3, \ 4$, by

where $D(\mathcal{A}_{h, 1}) = D(\mathcal{A}_{h, 2}) = D(\mathcal{A}_{h, 3}) = D(\mathcal{A}_{h, 4}) = {0_{L^1(\mathbb{R}+)}}\times W^{1, 1}(\mathbb{R}_+)$.

We can choose $\lambda\in\mathbb{C}$ such that $\mathit {Re} (\lambda)\geq-\mu_{0}$. Therefore, for $\lambda\in (r(\mathcal{A}_{h, 1})\cap r(\mathcal{A}_{h, 2})\cap r(\mathcal{A}_{h, 3})\cap r(\mathcal{A}_{h, 4}))$, and for any $(\psi_1, \psi_2)^T\in X_+^2$, the following explicit expression for the resolvent is available:

if and only if

Define $\mathcal{A}:D(\mathcal{A})\subset X\to X$ be a linear operator as

where, $D(\mathcal{A}) = D(\mathcal{A}_{h, 1})\times D(\mathcal{A}_{h, 2})\times D(\mathcal{A}_{h, 3})\times D(\mathcal{A}_{h, 4})\times D(\mathcal{A}_{v, 1})$. From the above discussion and Theorem 3.2 in ^[18], the linear operator $\mathcal{A}$ is a Hille–Yosida operator and the infinitesimal generator of $C_0$-semigroup. Further, let

and $\mathbb{X}_0 = \left(\{0_\mathbb{R}\}\times L^1(\mathbb{R}_+)\right)^5$, a nonlinear operator $\mathcal{F}:\mathbb{X}_0\to \mathbb{X}_0$ as

where, $\tilde{\lambda}_1(t, a) = \rho z_2(a)\int_{0}^{\infty}\beta_1(b)I_{v}(t, b)\mathrm{d}b$, $\tilde{\lambda}_2(t, a) = z_2(a)\int_{0}^{\infty}\beta_2(a)(\alpha A_{h}(t, a)+I_{h}(t, a))\mathrm{d}a$. Then, model (2.3) can be reformulated as the abstract Cauchy problem

From Theorem 3.5 in ^[15], or Theorem 3.2 in ^[18], the following result on the existence and uniqueness of a solution for the Cauchy problem (2.4) is valid.

Theorem 1. The problem (2.4) admits a unique global classical solution on $\mathbb{X}_0$; That is, model (2.3) has a unique global positive solution for the positive initial value.

Adding up the first four equations in model (2.3) gives the overall equation

with $N_{h}(t, 0) = \int_{0}^{\infty} b_{h}(a) N_{h}(t, a) \mathrm{d} a$, $N_{h}(0, a) = N_{h0}(a) = S_{h0}(a)+A_{h0}(a)+I_{h0}(a)+R_{h0}(a)$. System (2.5) is resolved following the characteristic curve $t-a = c$ (c is a constant)

To ensure the existence of a steady state, it is assumed that the population's net reproduction rate is identical to 1, i.e., $\int_{0}^{\infty} b_{h}(a)\mathrm{e}^{-\int_{0}^{a} \mu_{h}(\tau) \mathrm{d} \tau } \mathrm{d} a = 1.$ Hence, the steady state of system (2.5) is $N^{\infty}_{h}(a) = N^{\infty}_{h}(0)\mathrm{e}^{-\int_{0}^{a} \mu_{h}(\tau) \mathrm{d} \tau }$, where $N^{\infty}_h(a) = S_0(a)+A_0(a)+I_0(a)+R_0(a)$. By a simple calculation, we obtain

To simplify the initial boundary value problem, model (2.3) is normalized by

Then, model (2.3) and the force of infections become

and

with the initial and boundary conditions $s_{h}(t, 0) = 1$, $e_{h}(t, 0) = i_{h}(t, 0) = r_{h}(t, 0) = 0$ and

where $s_{h}(t, a)+e_{h}(t, a)+i_{h}(t, a)+r_{h}(t, a) = 1$.

3.   Stability of the disease-free steady state

There exists a disease-free steady state $\mathcal{E}^0 = (1, 0, 0, 0, 0)$ for model (2.6). Denote $\lambda_{1}(t, a)+\lambda_{2}(t, a) = z_2(a)W(t)$, where

To study the local stability of the disease-free steady state, it suffices to calculate the linearized system of model (2.6) at $\mathcal{E}^0$. We introduce the variable transformations $s_{h}(t, a) = \bar{s}_h(t, a)+1, e_{h}(t, a) = \bar{e}_h(t, a), i_{h}(t, a) = \bar{i}_h(t, a), r_{h}(t, a) = \bar{r}_h(t, a), i_{v}(t, b) = \bar{i}_v(t, b),$ and get a linearized system as

with the initial and boundary conditions

where

Assume that system (3.1) has the solution of exponential form $\bar{s}_h(t, a) = \bar{s}_h(a)\mathrm{e}^{\lambda t}$, $\bar{e}_h(t, a) = \bar{e}_h(a)\mathrm{e}^{\lambda t}$, $\bar{i}_h(t, a) = \bar{i}_h(a)\mathrm{e}^{\lambda t}$, $\bar{r}_h(t, a) = \bar{r}_h(a)\mathrm{e}^{\lambda t}$ and $\bar{i}_v(t, b) = \bar{i}_v(b)\mathrm{e}^{\lambda t}$, from which the following equation is obtained:

where

Solving the second, third and fifth equations of system (3.2) yields

Substituting $\bar{e}_h(a)$ in (3.4) into $\bar{i}_h(a)$ and $\bar{i}_v(b)$ yields

Substituting $\bar{e}_h(a)$ and (3.5) into (3.3), one gets

Dividing both the left and right sides of the above equation by $W_0$ ($W_0 \neq 0$), it follows that

The basic reproduction number of model (2.6) is defined as $\mathcal{R}_0 = :F(0)$, or expressly as

For the stability of the disease-free steady state $\mathcal{E}^0$, we obtain the following result.

Theorem 2. The disease-free steady state $\mathcal{E}^0$ of model (2.6) is locally asymptotically stable if $\mathcal{R}_0 < 1$ and unstable if $\mathcal{R}_0 > 1$.

Proof. According to (3.6), it is easy to see that $\mathrm{d}F(\lambda)/\mathrm{d}\lambda < 0$, and $\lim_{\lambda \rightarrow +\infty} F(\lambda) = 0$, $\lim_{\lambda \rightarrow -\infty} F(\lambda) = +\infty$. Consequently, the characteristic equation (3.6) has a unique real root. If $\mathcal{R}_0 < 1$, i.e., $F(0) < 1$, then $F(\lambda) = 1$ has a unique negative real root $\lambda^*$. Now, we claim that all other roots have negative real parts. In fact, if there is a root $\lambda = x+\mathrm{i}y$ with $x\geq 0$ such that $F(\lambda) = 1$. Then

Since $F(\lambda)$ is a monotonically decreasing function with respect to $\lambda$, it follows that $\mathit{Re} (\lambda) = x \leq \lambda^*$. This implies that all complex roots of Eq (3.6) have negative real parts. Hence, the disease-free steady state $\mathcal{E}^0$ of model (2.6) is locally asymptotically stable if $\mathcal{R}_0 < 1$. Conversely, if $\mathcal{R}_0 > 1$, indicating $F(0) > 1$, Eq (3.6) possesses a unique positive real root. In such a scenario, $\mathcal{E}^0$ is unstable. This concludes the proof. □

Theorem 3. The disease-free steady state $\mathcal{E}^0$ of model (2.6) is globally asymptotically stable if $\mathcal{R}_0 < 1$.

Proof. Define $g(t, a) = (\lambda_{1}(t, a)+\lambda_{2}(t, a)) s_h(t, a)$; It follows from the inequality $s_h(t, a) \leq 1$ that

By integrating model (2.6) over the characteristic line defined as $t-a = constant$, we obtain

When $a < t$ and $b < t$ are as specified in Eq (3.8), we obtain

Using Eqs (3.8) and (3.10), one obtains

Define $G(a) = \limsup_{t \rightarrow \infty} g(t, a)$ and, by applying Fatou's lemma ^[28] to the left and right ends of inequality (3.11) when $t \rightarrow \infty$, observe that

Assigning the constant $V$ as

Then $G(a) \leq z_2(a)V$, therefore

By the inequality, if $\mathcal{R}_0 < 1$, we conclude that $V = 0$, implying $G(a) = 0$, or equivalently, $\limsup_{t \rightarrow \infty}g(t, a) = 0$. Considering the expression in (3.8), one receives

Note that $s_h(t, a) = 1-e_h(t, a)-i_h(t, a)-r_h(t, a)$; Taking the limit, we obtain $\lim_{t \rightarrow \infty} s_h(t, a) = 1,$ as $t \rightarrow \infty$. Moreover, according to Theorem 1, if $\mathcal{R}_0 < 1$, the disease-free steady state $\mathcal{E}^0$ of model (2.6) is globally asymptotically stable. This concludes the proof. □

4.   Existence and stability analysis of the endemic steady state

Insights regarding the existence and stability of the endemic steady state $\mathcal{E}^*(s^*_h(a), e^*_h(a), i^*_h(a),$ $r^*_h(a), i^*_v(b))$ within model (2.6) can be derived as follows:

Theorem 4. If $\mathcal{R}_0 > 1$, there exists a unique endemic steady state $\mathcal{E}^*$ for model (2.6).

Proof. The endemic steady state $\mathcal{E}^*(s^*_h(a), e^*_h(a), i^*_h(a), r^*_h(a), i^*_v(b))$ of model (2.6) satisfies the following equations

where

Solve the first to fifth equations of (4.1)

Substitution $s^*_h(a)$ into $e^*_h(a)$ and $i^*_h(a)$ yield

Substituting (4.4) for $i^*_v(b)$ in (4.3) gives that

Substituting (4.4) and (4.5) into (4.2), and simplifying the equation yields

Letting $W^* = 0$ gives

Note that model (2.6) exhibits a unique endemic steady state $\mathcal{E}^*$, which corresponds to the existence of a unique positive root $W^*$ satisfying $H(W^*) = 1$. Considering the equation $s^*_h(a)+e^*_h(a)+i^*_h(a)+r^*_h(a) = 1$, where $s^*_h(a) > 0$ and $0 < \alpha < 1$, it implies that $\alpha e^*_h(a)+i^*_h(a) < 1$. Consequently, for any $W^* > 0$, there exists a

Considering the aforementioned inequality, there is $H(W^*) < 1$ when $W^* = \rho\int_{0}^{\infty}\beta_{1}(a)\mathrm{d} a+\int_{0}^{\infty}\beta_{2}(b)\mathrm{d} b$. Moreover, $H(W^*)$ behaves as a continuous and decreasing function of $W^*$, with $H(0) > 1$ for $\mathcal{R}_0 > 1$. Hence, equation (4.6) exists a unique real root on $(0, \rho\int_{0}^{\infty}\beta_{1}(a)\mathrm{d} a+\int_{0}^{\infty}\beta_{2}(b)\mathrm{d} b)$. This implies the existence of a unique $\mathcal{E}^*$ in model (2.6) when $\mathcal{R}_0 > 1$. The proof is concluded. □

The local asymptotic stability of the endemic steady state $\mathcal{E}^*$ is analysed below. If we consider $\tilde{s}_h(t, a)$, $\tilde{e}_h(t, a)$, $\tilde{i}_h(t, a)$, $\tilde{r}_h(t, a)$ and $\tilde{i}_v(t, b)$ as disturbances at the state $\mathcal{E}^*$, the linear system of model (2.6) at $\mathcal{E}^*$ takes the form

and

where

Moreover, assume that system (4.7) has an exponential solution of the form $\tilde{s}_h(t, a) = \tilde{s}_h(a)\mathrm{e}^{\omega t}$, $\tilde{e}_h(t, a) = \tilde{e}_h(a)\mathrm{e}^{\omega t}$, $\tilde{i}_h(t, a) = \tilde{i}_h(a)\mathrm{e}^{\omega t}$, $\tilde{r}_h(t, a) = \tilde{r}_h(a)\mathrm{e}^{\omega t}$ and $\tilde{i}_v(t, b) = \tilde{i}_v(b)\mathrm{e}^{\omega t}$. Thereby, the functions $\tilde{s}_h(a)$, $\tilde{e}_h(a)$, $\tilde{i}_h(a)$, $\tilde{r}_h(a)$ and $\tilde{i}_v(b)$ satisfy the following equations

where

It is important to acknowledge that the functions $\tilde{s}_h(a)$, $\tilde{e}_h(a)$, $\tilde{i}_h(a)$, $\tilde{r}_h(a)$, and $\tilde{i}_v(b)$ may have positive or negative values. Assuming $\hat {s}_h = {\tilde{s}_h}/{\tilde W}$, $\hat {e}_h = {\tilde {e}_h}/{\tilde {W}}$, $\hat {i}_h = {\tilde {i}_h}/{\tilde W}$, $\hat {r}_h = {\tilde {r}_h}/{\tilde W}$ and $\hat {i}_v = {\tilde {i}_v}/{\tilde W}$, then system (4.8) becomes

Based on the expression for $\tilde W$, it follows that

Solving the first to fifth equations of (4.9) yields

Substituting $\hat {e}_h(a)$ into $\hat {i}_h(a)$ and $\hat {i}_v(b)$ yields

Substituting (4.11), (4.12) into (4.10) yields

where

Proposition 1. Assume that $\varphi(a, \tau) > 0$ and $\psi(a, \tau) > 0$; For $\ 0\leq \eta \leq s \leq \tau \leq \xi \leq a$, then, $T(\omega)$ is a decreasing function of $\omega$ satisfies $\lim_{\omega \rightarrow +\infty} T(\omega) = 0$ and $\lim_{\omega \rightarrow -\infty} T(\omega) = +\infty$ and $T(0) < 1$.

Proof. By utilizing the expression in (4.13), ensuring that $\varphi(a, \eta) > 0$ and $\psi(a, \eta) > 0$ enables us to derive fundamental properties of $T(\omega)$, obtain $T(\omega)\geq 0$, $T'(\omega) < 0$, $\lim_{\omega \rightarrow +\infty} T(\omega) = 0$, $\lim_{\omega \rightarrow -\infty} T(\omega) = +\infty$. From Eq (4.13), letting $\omega = 0$ has

where

The proof is completed. □

The subsequent result pertains to the local asymptotic stability of the endemic steady state.

Theorem 5. Under $\varphi(a, \eta) > 0$ and $\psi(a, \eta) > 0$ for $\ 0\leq \eta \leq s \leq \tau \leq \xi \leq a$, then the endemic steady state of model (2.6) is locally asymptotically stable if $\mathcal{R}_0 > 1$.

Proof. Similar to the proof of Theorem 2, the above discussion and the Proposition 1, if $\mathcal{R}_0 > 1$, $\varphi(a, \tau) > 0$ and $\psi(a, \tau) > 0$ hold for all $0\leq \eta \leq s \leq \tau \leq \xi \leq a$, then there exists a unique negative real root $T(\omega) = 1$, and all complex roots possess negative real parts. Consequently, the endemic steady state of model (2.6) achieves local asymptotic stability if $\mathcal{R}_0 > 1$. This finishes the proof. □

5.   Uniform persistence

In this subsection, we consider the uniform persistence of model (2.6) using the continuation theory of an infinite dimensional dynamical system.

Lemma 1. Given that $u(t, \cdot; u_0) = (s_h(t, \cdot), e_h(t, \cdot), i_h(t, \cdot), r_h(t, \cdot), i_v(t, \cdot))$ represents the solution of model (2.6) with the initial condition $u_0 = (s_{h0}(a), e_{h0}(a), i_{h0}(a), r_{h0}(a), i_{v0}(b))\in X$, then

(i) there is a positive constant $\varepsilon_1 > 0$ such that $\liminf_{t\rightarrow \infty}\parallel s_h(t, \cdot)\parallel_{L^1} \geq \varepsilon_1$ for $s_{h0}(a) > 0$;

(ii) if there exists $t^*\geq 0$ such that $e_h(t^*, \cdot) > 0$, or $i_h(t^*, \cdot) > 0$, or $i_v(t^*, \cdot) > 0$, then $\varphi(t, \cdot) > 0$ for $t > t^*$, where, $\varphi(t, \cdot) = e_h(t, \cdot)$, $i_h(t, \cdot)$ and $i_v(t, \cdot)$.

Proof. Proof $(i)$, by $s_h(t, a)+e_h(t, a)+i_h(t, a)+r_h(t, a) = 1$ and $0 < \alpha < 1$, it follows that $\alpha e_h(t, a)+i_h(t, a) < 1$. Subsequently, based on model (2.6), it can be inferred that

where $\lambda(a) = z_2(a)\left[\rho\int_{0}^{\infty}\beta_{1}(b) \mathrm{d} b +\int_{0}^{\infty}\beta_{2}(a)\mathrm{d} a\right]$. Therefore, it is possible to obtain

Let's denote $\breve{s}^*_h(a) = \mathrm{e}^{- \int_0^a \lambda(s)\mathrm{d}s}$ and applying the comparison principle, one can conclude that $\liminf_{t\rightarrow \infty}s_h(t, a)\geq \breve{s}^*_h(a)$. Consequently, there is a $\varepsilon_1 > 0$ such that $\liminf_{t\rightarrow \infty}\parallel s_h(t, \cdot)\parallel_{L^1} \geq \varepsilon_1$.

Now, we turn to the (ii). If there exits a $t^* > 0$ such that $e_h(t^*, \cdot) > 0$, then one has $e_h(t, a) > 0$ for $t\geqslant t^*$ due to the expression of equation (3.7). Further, form the expressions of Eqs (3.8) and (3.9), one has $r_h(t, a) > 0$ and $i_v(t, b) > 0$ for $t\geqslant t^*$ due to the facts $i_h(t^*, a) > 0$ and $s_h(t^*, a) > 0$. Similarly, for $r_h(t^*, \cdot) > 0$ or $i_v(t^*, \cdot) > 0$, we also can get $\varphi(t, \cdot) > 0$ for $t\geqslant t^*$, here, $\varphi(t, \cdot) = e_h(t, \cdot)$, $i_h(t, \cdot)$ and $i_v(t, \cdot)$. This completes the proof of conclusion (ii). □

Theorem 6. If $\mathcal{R}_0 > 1$, for initial value $u_0 = (s_{h0}(\cdot), e_{h0}(\cdot), i_{h0}(\cdot), r_{h0}(\cdot), i_{v0}(\cdot))\in X$ with $e_{h0}(\cdot)+i_{h0}(\cdot)+i_{v0}(\cdot) > 0$, then there exists a constant $\varepsilon_2 > 0$ such that the solution $u(t, \cdot; u_0) = (s_h(t, \cdot), e_h(t, \cdot), i_h(t,$ $\cdot)$, $r_h(t, \cdot), i_v(t, \cdot))$ satisfies

Proof. Define

From Lemma 1, we understand that $Y$ serves as a positive invariant set for model (2.6) concerning the solution semi-flow $u(t)$, denote $N_\partial = \left\{u_0\in \partial{Y}:u(t, u_0)\in \partial{Y}, t\geqslant0\right\}$, where $\omega(u_0)$ denotes the omega limit set of $u(t, u_0)$. Now, one claims that $\{\mathcal{E}^0\} = \bigcup_{u_0 \in N_\partial} \omega(u_0)$.

It is evident that $u(t, \mathcal{E}^0) = \mathcal{E}^0$ holds for any $t\geq 0$, implying $\{\mathcal{E}^0\}\subset \bigcup_{u_0 \in N_\partial} \omega(u_0)$. On the contrary, for any $u_0 \in N_\partial$, we have $\int_{0}^{\underline{a}} e_h(t, a) \mathrm{d}a = 0$ or $\int_{0}^{\underline{a}} i_h(t, a) \mathrm{d}a = 0$ or $\int_{0}^{\underline{b}} i_v(t, b) \mathrm{d}b = 0$. Let's start by considering the scenario where $\int_{0}^{\underline{a}} e_h(t, a) \mathrm{d}a = 0$, implying $e_h(t, a)\equiv 0$. According to model (2.6), this leads to

Thus, the above equation holds if and only if $i_h(t, a) = i_v(t, b) = 0$. So there is $\left(\frac{\partial }{\partial t}+\frac{\partial }{\partial a}\right)s_h(t, a) = 0$ with $s_h(t, 0) = 1$. The $\lim_{t\rightarrow \infty}s_h(t, a) = 1$ by $s_h(t, a)+e_h(t, a)+i_h(t, a)+r_h(t, a) = 1$, and so $\omega(u_0) = \{\mathcal{E}^0\}$. Similar to the case of $\int_{0}^{\underline{a}} e_h(t, a) \mathrm{d}a = 0$, one also has $\omega(u_0) = \{\mathcal{E}^0\}$ when $\int_{0}^{\underline{a}} i_h(t, a) \mathrm{d}a = 0$ or $\int_{0}^{\underline{b}} i_v(t, b) \mathrm{d}b = 0$, thereby $\bigcup_{u_0 \in N_\partial} \omega(u_0) \subset \{\mathcal{E}^0\}$. In conclusion, it follows that $\{\mathcal{E}^0\} = \bigcup_{u_0 \in N_\partial} \omega(u_0)$. Therefore, all solutions of model (2.6) converge to $\{\mathcal{E}^0\}$ on $\partial{Y}$ as $t\rightarrow \infty$.

Further, it is shown that the uniformly weakly repulsive of $\{\mathcal{E}^0\}$, i.e., there exists a constant $\delta > 0$ such that $\limsup_{t\rightarrow \infty} \parallel u(t, u_0)-\mathcal{E}^0\parallel_X \geq\delta$ for any initial values $u_0\in X$. Using the counter-argument. Suppose that otherwise, that is, there exists an initial value $u_1\in X$ such that $\limsup_{t\rightarrow \infty} \parallel u(t, u_1)-\mathcal{E}^0\parallel_X < \delta$. Thus, there exists $\bar T > 0$ such that, for any $t > \bar T$,

It is clear that from model (2.6)

Get the auxiliary system as

with the initial and boundary conditions

where

The following exponential form of solution is obtained for system (5.1) as follows: $\underline{e_h}(t, a) = \underline{e_h}(a)\mathrm{e}^{pt}$, $\underline{i_h}(t, a) = \underline{i_h}(a)\mathrm{e}^{pt}$, $\underline{r_h}(t, a) = \underline{r_h}(a)\mathrm{e}^{pt}$ and $\underline{i_v}(t, b) = \underline{i_v}(b)\mathrm{e}^{pt}$, one has $\underline{\lambda_{1}}(t, a)+\underline{\lambda_{2}}(t, a) = z_2(a)\mathrm{e}^{p t} \Lambda,$ where

Consequently,

Solving the equations of system (5.2) yields

Substituting $\underline{e_h}(a)$ in (5.3) into $\underline{i_h}(a)$ and $\underline{i_v}(b)$ yields

Substituting (5.3) and (5.4) into the expression for $\Lambda$ yields

It is obvious from the expression for $D(p)$ that $D(p)$ is a monotonically decreasing function with respect to $p$, $\lim_{p\rightarrow \infty}D(p) = 0$, and

Taking a sufficiently small $\delta$ so that

It appears that when $\mathcal{R}_0 > 1$, there exists a unique positive root of Eq (5.5), meaning the solution $(\underline{e_h}(t, \cdot)$, $\underline{i_h}(t, \cdot)$, $\underline{r_h}(t, \cdot)$, $\underline{i_v}(t, \cdot))$ of system (5.1) becomes unbounded as $t > \bar T$. Consequently, the solution $(s_h(t, \cdot), e_h(t, \cdot)$, $i_h(t, \cdot)$, $r_h(t, \cdot)$, $i_v(t, \cdot))$ of model (2.6) also becomes unbounded for $t > \bar T$, which contradicts the boundedness of $u(t, \cdot; u_0)$. Hence, conclude that $\{\mathcal{E}^0\}$ is uniformly weakly repulsive.

To summarize, $\{\mathcal{E}^0\}$ is an isolated invariant set on $X$, and $W^s(\mathcal{E}^0)\bigcap Y = \emptyset$, where $W^s(\mathcal{E}^0)$ is a stable subset of $\mathcal{E}^0$. Moreover, there is no closed loop from $\mathcal{E}^0$ to $\mathcal{E}^0$ on $\partial Y$. By applying persistence theory ^[29], results in the uniform persistence of model (2.6). This completes the proof. □

6.   Optimal control problem

For all vector-borne infectious diseases, vaccination strategies for susceptible populations (denoted as $u_1(t, a)$) and elimination strategies for infectious vectors (denoted as $u_2(t)$) can be proposed. Therefore, by introducing the control variables $u_1$ and $u_2$ model (2.3) obtains

where

with the initial and boundary conditions

Considering practical scenarios, replace the upper limit of the age-related integral with a finite value $a^+ > 0$. the control set as follows

therefore, the objective function of our optimal control problem is

where $\tau_i$ and $B_i$ are defined as positive coefficients, serving to adjust the significance attributed to the state and control variables, respectively $(i = 1, \ 2, \ 3;\ j = 1, \ 2)$. The existence of optimal control is proofed as follows.

Theorem 7. There exist optimal control variables $u_1^*(t, a)$, $u_2^*(t, b)\in U$ such that

satisfies the initial and boundary conditions of system (6.1).

The demonstration of Theorem 7 closely resembles that of Theorem 3.1 in ^[25], hence, it will be omitted here.

Select an additional control $u_1^{\epsilon}(t, a) = u_1(t, a)+\epsilon l_1(t, a)$ and $u_2^{\epsilon}(t) = u_2(t)+\epsilon l_2(t)$, where $l_1(t, a)$ and $l_2(t)$ are variation functions and $\epsilon\in(0, 1)$, then

Therefore, the equation of state variables corresponding to the new control variables $u_i^{\epsilon}$ (i = 1, 2) are

where

with the boundary conditions

and the initial conditions

where $N^{\epsilon}_h(t, a) = S^{\epsilon}_h(t, a)+A^{\epsilon}_h(t, a)+I^{\epsilon}_h(t, a)+R^{\epsilon}_h(t, a)$, we find the following difference quotient

Where $\check{S}_h(t, a)$, $\check{A}_h(t, a)$, $\check{I}_h(t, a)$, $\check{R}_h(t, a)$ and $\check{I}_v(t, b)$ comply with the following system

with the boundary conditions

and the initial conditions

where $\check{N}_h(t, a) = \check{S}_h(t, a)+\check{A}_h(t, a)+\check{I}_h(t, a)+\check{R}_h(t, a)$. In order to find the adjoint equations, we consider the first equation of the system (6.3) as

with the conditions $\check{S}_h(0, a) = 0$, $\check{S}_h(t, a^+) = 0$, $\Lambda^*_1(T, 0) = 0$ and $\langle\langle f, g\rangle\rangle_1 = \int_0^T\int_0^{a^+}fg \mathrm{d}a\mathrm{d}t$. Similarly, the rest of the equations for system (6.3) can become

under the initial conditions $\check{A}_h(0, a) = 0$, $\check{A}_h(t, a^+) = 0$, $\Lambda^*_2(T, 0) = 0$.

under the initial conditions $\check{I}_h(0, a) = 0$, $\check{I}_h(t, A^+) = 0$, $\Lambda^*_3(T, 0) = 0$.

under the initial conditions $\check{R}_h(0, a) = 0$, $\check{R}_h(t, a^+) = 0$, $\Lambda^*_4(T, 0) = 0$.

under the initial conditions $\check{I}_v(0, b) = 0$, $\check{I}_v(t, b^+) = 0$, $\Lambda^*_5(T, 0) = 0$ and $\langle\langle f, g\rangle\rangle_2 = \int_0^T\int_0^{b^+}fg\mathrm{d}b\mathrm{d}t$.

Next the Lagrangian $\mathcal{L}$ function is defined as

By solving $\frac{\partial\mathcal{L}}{\partial\check{S}_h} = 0$, $\frac{\partial\mathcal{L}}{\partial\check{A}_h} = 0$, $\frac{\partial\mathcal{L}}{\partial\check{I}_h} = 0$, $\frac{\partial\mathcal{L}}{\partial\check{R}_h} = 0$, $\frac{\partial\mathcal{L}}{\partial\tilde{I}_v} = 0$, combined with Eqs (6.4)–(6.8), then we obtain the adjoint equations

with the transversality conditions

with the boundary conditions

Theorem 8. If $u_1^*$, $u_2^*\in U$ represent optimal controls that minimize the objective function $\mathcal{J}(u_1, u_2)$, and $(S^*_h(t, a), A^*_h(t, a), I^*_h(t, a), R^*_h(t, a), I^*_v(t, b))$ and $(\Lambda^*_1(t, a), \Lambda^*_2(t, a), \Lambda^*_3(t, a)$, $\Lambda^*_4(t, a), \Lambda^*_5(t, b))$ denote the corresponding state and adjoint variables, respectively, then

Proof. The Gateaux derivative of $\mathcal{J}(u_1, u_2)$ is

Using the optimal values of the state variables, the above inequality can be expressed as

Hence, the optimal control variable can be characterized as

where $u_i^*\in L^{\infty}(\Omega)$ and $0\leq u_i^*\leq h_i(i = 1, 2)$. This finishes the proof. □

Subsequently, we confirm the existence of a unique optimal control strategy. Given the complexity of finding control sequences and associated states that converge to the optimal controls and states in the partial differential equation model with age structure, we resort to obtaining minimizing sequences of approximate functions using Ekeland's Principle ^[30] from ^[31]. Let us assume the existence of a pair of control variables that minimize the following objective function:

Theorem 9. If $(u_1^{\epsilon}, u_2^{\epsilon})$ is the minimization sequence of $\mathcal{J}_{\epsilon}(u_1, u_2)$, then

where the function $\theta_i^{\epsilon} \in L^{\infty}(\Omega)$ such that $|\theta_i^{\epsilon}|\leq1(i = 1, 2)$ for all $(t, a)\in \Omega$.

The proof of this theorem is similar to Theorem 8 and is omitted here.

Theorem 10. There exists a unique optimal control $(u_1^*, u_2^*)$ to minimize objective function $\mathcal{J}(u_1, u_2)$ if $\frac{T}{B_1}$ and $\frac{T}{B_2}$ are sufficiently small.

Proof. Define two functions by

basing on ^[31], taking into account both sets of control variables $(u_1, u_2)$ and $(\hat{u}_1, \hat{u}_2)$, along with the Lipschitz properties of the state and adjoint variables, we can derive the following conclusions:

given the $L^{\infty}$ bounds of the state and its adjoint solution. Here, $K_1$ and $K_2$ are ascertained, correlating with these bounds and the Lipschitz constants. Provided that $\frac{T}{B_1}$ and $\frac{T}{B_2}$ are sufficiently small, then

This suggests that $(u_1^{\epsilon}, u_2^{\epsilon})$ converges to $(u_1^*, u_2^*)$. Applying Ekeland's principle, we can conclude that

The proof is completed. □

7.   Numerical simulations

To better interpret the theoretical results and to analyze the impact of age structure and multiple transmission routes on the spread of the disease, some numerical simulations are conducted. For this purpose, the finite difference method is used to discretize the extended eigenline of model (2.3) as

where $C_1 = \int_{0}^{\infty}\beta_1(b)I_{v}(t, b)\mathrm{d}b$, $C_2 = \int_{0}^{\infty}\beta_2(a)(\alpha A_{h}(t, a)+I_{h}(t, a))\mathrm{d}a$, and use the complex trapezoidal formula for linear approximations. Then fix some basic parameters of model (2.3) as $\alpha = 0.45$, $q = 0.6$, $k(a) = 1.5\times 10^{-4}(1+\sin(0.1\pi a))$, $\gamma_1(a) = 0.38(1+\sin(0.1\pi a))$, $\gamma_2(a) = 0.16(1+\sin(0.1\pi a))$, $\beta_3(a) = 4\times 10^{-5}(1+\sin(0.1\pi a))$, $\mu_h(a) = \frac{0.01}{80-a}$.

First, the effects of the probabilities of human-to-human (that is, $\beta_2(a)$) and vector-to-human (that is, $z_1(a)$) on the basic reproduction number $\mathcal{R}_0$ are discussed. According to the expression for the reproduction number $\mathcal{R}_0$, $\beta_2(a)$ and $z_2(a)$ are positively correlated with $\mathcal{R}_0$, which is shown in Figure 1(a). In addition, from Figure 1(b), it also reveals that when $\beta_2(a) < 0.1$, $z_2(a) = 0.1(1+\sin(0.1\pi a))$, then $\mathcal{R}_0 < 1$. Therefore, appropriate use of some preventive measures (wearing long clothes and long sleeves, vaccination, going to fewer places where people congregate, etc.) can lead to a reduction in $\mathcal{R}_0$, thus controlling the spread of vector-borne diseases.

Now, we choose $\beta_2(a) = 0.05(1+ \sin(0.1\pi a))$, $z_2(a) = 0.1(1+\sin(0.1\pi a))$, $\beta_1(b) = 7.5\times 10^{-4}(1+\sin(0.1\pi b))$, and $\beta_2(a) = 1.1\times 10^{-4}(1+ \sin(0.1\pi a))$, the basic reproduction number $\mathcal{R}_{0}\approx 0.4893 < 1$ by direct calculation. From Theorem 3, we know that the disease-free steady state of model (2.6) is globally asymptotically stable. In this scenario, the density distributions of the asymptomatic infected individuals and infected vectors are shown in Figure 2(a) and (b), respectively. That is, the number of infected individuals and vectors of infection in all age groups gradually approaches $0$ as time $t$ increases. In addition, Figure 2(c) and (d) show that although the disease is ultimately extinct, the distribution of infected individuals is still age-heterogeneous. Therefore, in the process of disease prevention and control, the limited medical resources can be more reasonably deployed by fully considering the age distribution characteristics of the infected.

However, if we change the transmission rates to $z_2(a) = 0.3(1+\sin(0.1\pi a))$, $\beta_1(b) = 0.0025(1+\sin(0.1\pi b))$ and $\beta_2(a) = 0.03(1+\sin(0.1\pi a))$, then $\mathcal{R}_{0}\approx 2.606 > 1$. According to Theorem 4, exists the unique endemic steady state for model (2.6). The plots in Figure 3(a) and (b) represent that the distributions of symptomatic infected individuals and infected vectors, respectively. Meanwhile, the numerical simulation results indicated that the period and size of the disease outbreaks varied among different age groups and are very closely related to the actual age of the individuals and the age of infection of the vectors. Therefore, the age structure factor plays a crucial role in the spread of the vector-borne disease. As can be seen from the plots, $I_h(t, a)$ and $I_v(t, b)$ tend to the endemic steady state as $t$ tends to infinity for different initial values. In the same manner, as time approaches infinity, the density distributions of $S_h(t, a)$, $A_h(t, a)$, and $R_h(t, a)$ converge towards their respective endemic steady state, indicating asymptotic stability. This implies that the disease is persistent.

Finally, we consider the effect of the control strategy on the distribution of disease. For this purpose we choose all control strengths to be constants (that is, $u_1(t, a) = u_1$ and $u_2(t) = u_2$ are constants) and focus on characterizing the age distribution of those who contract the disease for the same control intensity. The plots in Figure 4(a) and (b) show the effect of the control strategy on symptomatic infected individuals and infected vectors, which shows that the total number of infected individuals can be significantly reduced by increasing vaccine coverage (or awareness of personal protection) in susceptible individuals and insecticide spraying of infected vectors. Further, the red and blue curves in Figure 4(c) and (d) indicate the age distribution of infected hosts and infected vectors, respectively, in the early stages of control (i.e., when $t$ is relatively small). Numerical simulations show that vaccination rates (or personal protective measures) provide better protection for younger populations, while mosquito eradication rapidly reduces the number of infected vectors. More fundamentally, when discussing the impact of control measures on disease transmission, as depicted in Figure 5(a), the peaks of infected hosts progressively diminished, and the total number of infected hosts significantly decreased with the escalation of control measure $u_1$ intensity from 0.1 to 0.3, and then to 0.9. Additionally, Figure 5(b) illustrates that the number of infected vectors also gradually decreased with the increasing intensity of control measure $u_2$, which was varied from 0.15 to 0.25 and ultimately to 0.95. Consequently, it is evident that the use of high-intensity control measures can lead to a significant reduction in both the number of infected individuals and vectors.

8.   Conclusions

Considering the prevalence of asymptomatic infections and the distinct differences in social activities among different age groups, in this paper, we construct a novel model that incorporates physiological age and factors such as multiple transmission routes and asymptomatic infections, setting our model apart from existing ones. Using linear approximation, the comparison principle, and Fatou's lemma, we derive the basic reproduction number, denoted as $\mathcal{R}_{0}$. We demonstrate that the disease-free steady state is globally asymptotically stable if $\mathcal{R}_{0} < 1$. Conversely, under specific conditions, a unique endemic steady state exists and is locally asymptotically stable when $\mathcal{R}_{0} > 1$. Furthermore, when $\mathcal{R}_{0} > 1$, the disease exhibits uniform persistence, which is a key focus of this paper.

In addition, we introduce two control measures, $u_1(t, a)$ and $u_2(t)$, representing vaccination or personal protection for susceptible individuals and insecticide spraying for vectors, respectively, aimed at disease control. This extension transforms our model into an optimal control problem involving partial differential equations with multiple age structures. We establish the existence and uniqueness of solutions to the optimal control problem and the associated equations for the control variables by employing the Gâteaux derivative. This approach is uncommon in the existing literature and represents a significant contribution of this paper. Finally, the theoretical findings of this paper are substantiated by numerical simulations, which confirm the impact of multiple transmission routes and age-structured factors on the spread of vector-borne disease.

Nevertheless, there remains a substantial territory that warrants further exploration. While we have delineated specific conditions for establishing the local asymptotic stability of the endemic steady state, our numerical simulations suggest that these conditions may be unnecessary. Therefore, it is essential to investigate methods for achieving local asymptotic stability of the endemic steady state without imposing additional constraints. Furthermore, to conduct a rigorous theoretical analysis, this study has excluded the influences of climate change variables—such as temperature, humidity, and rainfall—as well as uncertainties regarding vector population size and behavior in the mathematical modeling. Consequently, developing vector-borne disease models that incorporate multifactorial dynamics and examining how various factors influence disease prevention and control strategies represent compelling areas for future research.

Author contributions

Huihui Liu: Formal analysis, Conceptualization, Writing-original draft, Writing–review and editing; Yaping Wang and Linfei Nie: Funding acquisition, Supervision, and Editing. All authors have accepted responsibility for the entire content of this manuscript and consented to its submission to the journal, reviewed all the results, and approved the final version of the manuscript.

Use of Generative-AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors would like to thank the anonymous reviewers and the editors for their valuable suggestions for the improvement of the paper. This research is partially supported by the Tianshan Talent Training Program (Grant No.2022TSYCCX0015), the National Natural Science Foundation of China (Grant No. 12361103), and the Scientific Research and Innovation Project of Outstanding Doctoral Students in Xinjiang University (Grant No. XJU2024BS044).

Conflict of interest

The authors declare that there is no conflict of interest regarding the publication of this paper.