Empirical likelihood for varying coefficient partially nonlinear model with missing responses

1.   Introduction

Human immunodeficiency virus (HIV) infection was regarded as a major global public health issue. The World Health Organization (WHO) has claimed more than 40.4 million (32.9–51.3 million) infection cases as of 13 July 2023, in which 0.63 million (0.48–0.88 million) individuals died from HIV-related causes during the period from January 2022 to December 2022 ^[1]. Moreover, China admitted 117.9 thousand HIV/AIDS infection cases with 30.0 thousand ($25.48\%$) deaths from January 2007 to December 2011 ^[2], in which 3266 HIV/AIDS infection cases were diagnosed with 725 deaths ($22.20\%$), among which, 2320 local HIV/AIDS infection cases with 530 deaths ($22.84\%$) and, 946 HIV/AIDS infection cases with 195 deaths ($20.61\%$) from other provinces and overseas were reported during the period from January 2004 to December 2011. Additionally, the symptom characteristics were well recorded and classified by the surveillance data. Infection cases with mild symptoms may have an influenza-like illness including fever, headache, rash, sore throat and other symptoms such as swollen lymph nodes, weight loss, diarrhea and a cough, as reported by the WHO. Without treatment, infection cases with mild symptoms can develop to include severe symptoms such as tuberculosis, cryptococcal meningitis, severe bacterial infections and cancers due to the weakness of the immune system ^[1]. Essentially, the severity of symptoms was characterized by CD4$^+$ T cell counts ^[3,4,5], and the starting and ending points of structured treatment interruptions were determined in ^[6]. The transmission mechanism of symptom-dependent HIV/AIDS was usually described by mathematical models.

Local governments and policy makers controlled and prevented infectious diseases by adopting mathematical models, which played irreplaceable roles and provided new perspectives when studying the complex dynamics of infectious diseases, as in ^[7]. Specifically, the compartment models were often governed to investigate the transmission mechanisms of infectious diseases, in which the compartment structures were determined by specific infectious diseases. In fact, various types of deterministic dynamical models were included such as ordinary differential equation models, difference equation models, delay-differential equation models, age-structured partial differential equations (PDE) models and diffusion models in ^{[7,8,9,10,11,12,13,14,15,16]}. Meanwhile, recent contributions also governed the stochastic differential equation models and the fractional differential equation models to investigate the basic reproduction number and dynamical properties in ^{[17,18,19,20,21,22,23,24,25]}. More precisely, HIV/AIDS incorporated with tuberculosis (TB) models were taken into account, where the effect of Antiretroviral Therapy (ART) treatment and both local and global stabilities of HIV-only model were studied in ^[19,20] by adopting the standard incidence rates, in which the transmission mechanism to the susceptible compartment in ^[19] presented a more complex incidence mechanism. After these two contributions, the parameter fluctuation was introduced into the HIV/AIDS model with the bilinear incidence rate in ^[21], and the systematic fluctuation was investigated in ^[22] based on the model in ^[20]; therein, the related survival analysis and dynamical properties were extensively studied. The research results of the stochastic HIV/AIDS models in ^{[21,22,25,26]} showed that the intensities of white noises diminished the scale and threshold of HIV/AIDS. The investigations of HIV/AIDS models in ^[23,24] reflected that protection awareness of the susceptible also diminished the number of HIV/AIDS infection cases.

In this paper, we proposed symptom-dependent HIV/AIDS models to describe the epidemiological characteristics and dynamical properties of HIV/AIDS in the Fujian Province. Based on the features of real surveillance data from the Fujian Provincial Center for Disease Control and Prevention (Fujian CDC), the threshold $R_0$ of the deterministic HIV/AIDS model was derived, and the local and global stabilities for the disease-free point and endemic equilibrium point were extensively investigated in Section 3. Meanwhile, we introduced environmental fluctuations into the deterministic HIV/AIDS model; the conditions for the stationary distribution and stochastic extinction were investigated in Section 4. Moreover, we conducted numerical simulations and made predictions on the scale of HIV/AIDS infections in the Fujian Province in Section 5. The results of this study provided new perspectives to control the scale of HIV infection for local governments and policy makers including, but not limited to, the Fujian Province.

2.   Model formulation

HIV/AIDS infection cases from surveillance data (2004–2011) in the Fujian CDC were recorded as either mild cases or severe cases. Based on the transmission mechanism and symptoms of HIV/AIDS, in this paper, we proposed symptom-oriented HIV/AIDS models in a local population. Here, the total population was separated into four the following mutually-exclusive compartments: $S$, susceptible with no HIV/AIDS infection; $E$, HIV infected with no clinical symptoms but were able to transmit the HIV virus to others (HIV virus lived in the hosts but did not produce clinical symptoms, the hosts with an HIV virus lacked awareness of going to hospital and getting checked); $I_m$, HIV infected with mild symptoms; and $I_s$, HIV infected with severe symptoms. Hence, the total population at time $t$ for a designated region was given by $N(t) = S(t)+E(t)+I_m(t)+I_s(t).$ Especially, for the key men who have sex with men (MSM) population with HIV/AIDS, the transmission mechanism of HIV/AIDS was usually described by the bilinear incidence between $S$ and $E$, as well as $S$ and $I_m$ in recent contributions such as equations-oriented descriptions in ^{[21,23,25,27]} and fluctuation-oriented descriptions in ^[28,29]. The aforementioned HIV/AIDS models with extensive discussions motivated us to govern the bilinear incidence for investigating symptom-dependent HIV/AIDS models in this paper:

where $\Lambda$ was the constant recruitment rate, $\beta_e$ depicted the effective contact coefficient between $S$ and $E$, $\beta_m$ was the effective contact coefficient between $S$ and $I_m$, $1/\gamma_e$ represented the average time from the date of being infected by the HIV virus to the date with mild symptoms, $1/\gamma_m$ was the average time required from the date of being detected with mild symptoms to the date with severe symptoms, $1/\delta$ was the average time from the date that they had severe symptoms to the date that they died, and $\mu$ denoted the nature death rate. All parameters were assumed to be positive by their biological meanings.

3.   Survival investigation of model (2.1)

Let $X = (S, E, I_m, I_s)^\text{ T}$. Adding up the four equations in (2.1), and using the standard comparison theorem, we obtained that $N(t)\leq \lambda/\mu$; furthermore, we considered the dynamical properties in the positive invariant sets $\Omega = \left\{X \in \mathbb R_+^4 | 0 \leq S+E+I_m+I_s \leq \Lambda/\mu\right\}$ and $\Omega_0 = \{X\in \Omega | E = I_m = I_s = 0\}$ in this paper. It was easy to check that $P_0 = \left(S_0, E_0, I_{m0}, I_{s0}\right)^\text{ T} = \left(\Lambda/\mu, 0, 0, 0\right)^\text{ T}$ was the disease-free equilibrium point. It was known that the basic reproduction number was an important threshold for describing the average number of new HIV/AIDS infection cases produced by one HIV infected individual, which was computed by the next generation matrix in ^[30]. Let $\mathscr F$ be the new infections in compartments $E, I_m, I_s$, and let $\mathscr V = \mathscr V^- - \mathscr V^+$, in which $\mathscr V^-$ was the transfer rate for the individuals moving out of three compartments, $\mathscr V^+$ was the transfer rate for the individuals entering into the compartments by all other means. Then, we wrote down the following:

The Jacobian matrices for $\mathscr F$ and $\mathscr V$ were computed respectively as follows:

Substituting the disease-free equilibrium point $P_0$ into matrix $FV^{-1}$ gave the following basic reproduction number:

where $\rho(\cdot)$ denoted the spectral radius, and $k_1 = \gamma_m+\mu, \quad k_2 = \delta+\mu, \quad k_3 = \gamma_e+\mu, \quad k_4 = \gamma_e \beta_m + k_1 \beta_e.$ $R_0$ reflected the average number that an HIV/AIDS infection case transmitted the HIV virus to the susceptible individuals without interventions.

Alternatively, the endemic equilibrium point $P^* = (S^*, E^*, I_m^*, I_s^*)^\text{ T}$ existed for model (2.1) when $R_0 > 1$, where

In fact, let $f\in\mathbb R$ denote the right-hand side of the second equation in (2.1) and substitute $S^*, I_m^*, I_s^*$ into the equation. Then, the following expression was obtained:

in which $A = k_3k_4 > 0, \; B = k_4 \Lambda-k_1k_3\mu = k_1k_3\mu(R_0-1) > 0$, $\hat f(E^*) = -AE^{*2}+BE^*$, $g(E^*) = k_4 E^*+k_1\mu.$ It was easy to check that $\hat f(E^*)$ is a quadratic function with a negative quadratic coefficient and $f(0) = 0$; therefore, when $R_0 > 1$, then $\hat f(E^*) = 0$ had a unique positive real root, which implied that $f(E^*) = 0$ also admitted a unique positive real root since $g(E^*) > 0$ held for $E^* \in \mathbb R^+$.

Furthermore, in this part, the locally and globally asymptotic stabilities for two equilibrium points, $P_0$ and $P^*$, to model (2.1) were concerned using the Routh-Hurwitz Theorem and LaSalle's Invariant Principle, respectively.

Theorem 3.1. ($i$). If $R_0 < 1$, then the disease-free equilibrium point $P_0$ was locally asymptotically stable in domain $\Omega$.

($ii$). If $R_0 > 1$, then $P_0$ was unstable, and the endemic equilibrium point $P^*$ was locally asymptotically stable in $\Omega \backslash \Omega_0$.

Proof. (ⅰ) The Jacobian matrix of model (2.1) was expressed as follows:

Substituting $P_0$ into $J(X)$, then the eigenvalues of the characteristic equation $\det (\lambda \mathcal{I}-J(P_0)) = 0$ were given by $\lambda_1 = -\mu < 0, \lambda_2 = -k_2 < 0$ and

Here,

If $\beta_e S_0 - \gamma_e + \gamma_m \geq 0$, then $\lambda_3 < -\gamma_m-\mu < 0;$ if $\beta_e S_0 - \gamma_e + \gamma_m < 0$, then we obtained that $\lambda_3 < k_3(R_0-1)-\frac{\gamma_e\beta_m S_0}{k_1} < 0$ by $R_0 < 1$. The product of $\lambda_3$ and $\lambda_4$ gave the following expression:

together with $\lambda_3 < 0$, which implied that $\lambda_4 < 0$. By the Routh-Hurwitz Theorem, the disease-free equilibrium point $P_0$ was local asymptotically stable.

(ⅱ). If $R_0 > 1$, then $\lambda_3\lambda_4 = -k_1k_3(R_0-1) < 0$. Thus, there existed one positive eigenvalue, and the disease-free equilibrium point $P_0$ was unstable.

Then, we substituted $P^*$ into $J(X)$ to investigate the local asymptotic stability of the endemic equilibrium point $P^*$ as follows:

Obviously, $\lambda_4 = -k_2 < 0$. We denoted $k_5 = k_3(k_4E^*+k_1\mu) = \Lambda k_4.$ Let the characteristic equation of $J(P^*)$ be $\lambda^3 + a_2\lambda^2 + a_1\lambda + a_0 = 0,$ where

Then,

By the Routh-Hurwitz Theorem, the endemic equilibrium point $P^*$ was locally asymptotically stable in $\Omega \backslash \Omega_0$. □

Theorem 3.2. ($i$). If $R_0 \leq 1$, then the disease-free equilibrium point $P_0$ was globally asymptotically stable in domain $\Omega$.

($ii$). If $R_0 > 1$, then the endemic equilibrium point $P^*$ was globally asymptotically stable in $\Omega \backslash \Omega_0$.

Proof. Noting that the compartment $I_s$ did not appear in the first three equations of model (2.1), we considered the following equivalent model:

with the following positive invariant sets:

Then, the dynamics of model (2.1) was the same as model (3.2).

$(i)$. We defined a $C^2$-function $V_1:\mathbb{R}_+^3\rightarrow \mathbb{R}_+$ by the following:

The time derivation of $V_1$ along the positive solution was given by the following:

Using relation $S_0 = \frac{\Lambda}{\mu}$, we derived the following:

Thus, $\dot V_1 \leq 0$ when $R_0 \leq 1$. The equality held if and only if $S = S_0, E = 0$, which corresponded to $\overline\Omega_1 = \left\{(S, E, I_m)\in \overline\Omega:S = \Lambda/\mu, E = 0\right\} \subset \overline\Omega.$ If $R_0 = 1$, then $\dot V_1 \leq 0$; the equality was valid if and only if $S = S_0$, which corresponded to $\overline\Omega_2 = \left\{(S, E, I_m)\in \overline\Omega:S = \Lambda/\mu\right\} \subset \overline\Omega.$ Recalling the boundedness of $\overline\Omega,$ the maximal compact invariant set of model (3.2) on $\overline\Omega_1$ and $\overline\Omega_2$ contained only one element, $P_0$. According to LaSalle's Invariant Principle, the disease-free equilibrium point $P_0$ was globally asymptotically stable.

$(ii)$. For convenience, we denoted $x = \frac{S}{S^*}, \; y = \frac{E}{E^*}, \; z = \frac{I_m}{I_m^*}.$ By multiplying $\frac{1}{S^*}, \frac{1}{E^*}$ and $\frac{1}{I_m^*}$ on both sides of model (3.2), together with the identities

model (3.2) was transformed into the following form:

We constructed a $C^2$-function $V_2:\mathbb{R}_+^3\rightarrow \mathbb{R}_+$ by

to verify the global stability of the endemic equilibrium point $\bar P^*$ for model (3.3). Differentiating $V_2$ with respect to $t$ along the positive solution of model (3.3), we obtained the following:

The following relationships $\Lambda - \beta_mS^*I_m^* > \Lambda - \beta_eS^*E^* - \beta_mS^*I_m^* - \mu S^* = 0, \beta_mS^*I_m^* > 0$ and $2-x-\frac{1}{x} \leq 0, 3-\frac{1}{x}-\frac{xz}{y}-\frac{y}{z} \leq 0$ yielded that $\dot V_2 \leq 0$. The equality was valid if and only if $x = 1$ and $y = z$, which implied that $S = S^*$ and $\frac{E}{E^*} = \frac{I_m}{I_m^*}$. Therefore, the maximal compact invariant set of model (3.2) on

was a singleton $\{P^*\}$. Hence, by LaSalle's Invariant Principle, the unique endemic equilibrium point $P^*$ was globally asymptotically stable for $R_0 > 1$. □

4.   Survival investigation of model (4.1)

Epidemic models with environmental fluctuations were profitable to provide additional perspectives, such as the survival analysis, as compared with their deterministic counterparts ^[31]. In general, environmental fluctuations contained linear fluctuations in ^[32,33] and nonlinear fluctuations in ^[34]. Some studies have shown that the evolving process of HIV was naturally subject to environmental fluctuations ^[21,22,25] such as policies, medical systems, climate and so on. Motivated by the previous contributions, in this paper, we assumed that $X = (S, E, I_m, I_s)^\text{ T}$ was a Markov process ^[35] and assumed that the environmental fluctuations were proportional to $S, E, I_m$, and $I_s$. Taking the environmental fluctuations into account, model (4.1) was described by the following form:

where $B_i(t)$ were four independent standard Brownian motions defined on a complete filtered probability space $(\Omega, \mathcal F, \{\mathcal F_t\}_{t \geq 0}, \mathbb P)$ with a filtration $\{\mathcal F_t\}_{t \geq 0}$, which was increasing and right continuous while $\mathcal F_0$ contained all $\mathbb P$-null sets ^[35] with the initial values $B_i(0) = 0$; $\sigma_i$ reflected the intensities of white noises for $i = 1, 2, 3, 4$.

Then, we showed that there existed a unique global positive solution to model (4.1) for any given initial values, which was described by the undermentioned Theorem 4.1.

Theorem 4.1. Model (4.1) had a unique global positive solution $X(t) \in \mathbb{R}_+^4$ with the initial value $X_0\in \mathbb{R}_+^4$ for any $t\geq 0$.

Proof. Since the coefficients of model (4.1) were locally Lipschitz continuous, there existed a unique local solution $X(t)$ on $t\in[0, \tau_e)$ for any initial value $X_0 \in \mathbb{R}_+^4$, where $\tau_e$ denoted the explosion time. In order to prove that $X(t)$ was global, we needed to verify that $\tau_e = \infty$ held almost surely. Let $n_0 > 1$ be sufficiently large such that each component of $X(0)$ stayed in $[{\frac{1}{n_0}}, n_0]$. Let the infimum of an empty set equals $\infty$. Obviously, $\left\{\tau_n\right\}_{n\geq n_0}$ was monotonically increasing as $n\rightarrow \infty$. Set $\tau_\infty = \lim_{n \to \infty}{\tau_n}$; then, we obtained $\tau_\infty\le\tau_e$ by the definition of stopping time. The proof was given by a contradiction. We assumed that there existed a pair of positive constants $T > 0$ and $\varepsilon\in(0, 1)$ such that the probability that $\tau_\infty\le T$ was larger than $\varepsilon$ for any $n\geq n_0$ (i.e., $\mathbb{P}\left\{\tau_n\le T\right\}\geq\varepsilon$). We defined a $C^2$-function $V_1:\mathbb{R}_+^4\rightarrow \mathbb{R}_+$ by the following:

Applying the Itô's formula, we obtained the following:

The remaining parts of the proof followed the approaches in ^[36,37,38] and we omitted them here.

Next, we provided two useful results, Lemmas 4.1 and 4.2, and the corresponding proofs were quite similar with Lemmas 2.1 and 2.2 in ^[39]; therefore, we omitted them here. We denoted $X_0 = (S(0), E(0), I_m(0), I_s(0))^\text{ T}.$

Lemma 4.1. Let $X(t)$ be a solution of model (4.1) initiated with $X_0\in \mathbb{R}_+^4$; then,

Lemma 4.2. Suppose that $\mu > (\sigma_1^2\vee \sigma_2^2\vee \sigma_3^2 \vee \sigma_4^2)/2$. Let $X(t)$ be a solution of model (4.1) initiated with $X_0\in \mathbb{R}_+^4$; then,

Theorem 4.2. For any initial value $X_0\in \mathbb{R}^4_+$, the solution of model (4.1) was stochastically ultimately bounded.

Proof. We defined $V_2(X) = S^\theta + E^\theta + I_m^\theta + I_s^\theta,$ where $\theta \in (0, 1)$. By using Itô's formula, we obtained $\mbox{d}V_2(X) = \mathcal{L}V_2(X)\mbox{d}t + G(X)\text{d}B(t),$ in which

and

We defined $F(X) := V_2(X)+\mathcal{L}V_2(X),$ which was bounded in $\mathbb{R}^4_+$ to model (4.1). Furthermore, there existed a constant $H_1$ such that $F(X)\leq H_1 < \infty$; hence, $\mathcal{L}V_2(X)\leq H_1-V_2(X).$ Again, by Itô's formula, for $e^tV_2(X)$, we obtained the following:

Integrating (4.2) and taking the expectation yielded the following:

Let $n\rightarrow \infty$, we obtained the following:

Noticing that

we obtained

which revealed that

Now, for any $\varepsilon > 0$, let $H = (2^\theta H_1/\varepsilon)^2$. By Chebyshev's Inequality, we obtained the following:

Let $\theta = 1/2$. Then, we derived

which was rewritten as follows:

The proof was complete. □

In the undermentioned Theorem 4.3, we established the sufficient conditions to guarantee the existence of a stationary distribution and the ergodicity of model (4.1). This main result showed that HIV exhibited a sustainable behavior in a long run. We denoted

which degenerated to $R_0$ when $\sigma_1 = \sigma_2 = \sigma_3 = 0$. Here, the expression of $R_0^s$ was independent of the fluctuation of $I_s$. In order to show the main result for Theorem 4.3, we used Lemma 3.1 of ^[25] to check the details.

Theorem 4.3. Model (4.1) admitted a unique stationary distribution, and it was ergodic when $R_0^s > 1$.

Proof. The positive definite diffusion matrix of model (4.1) was $A\left(X\right) = \text{diag} \left\{\sigma_1^2S^2, \sigma_2^2E^2{, \sigma}_3^2I_m^2{, \sigma}_4^2I_s^2\right\};$ therefore, condition (ⅰ) of Lemma 3.1 of ^[25] was satisfied. To prove condition (ⅱ), we constructed a $C^2$-function $J\left(X\right) := MV_3+V_4+V_5+V_6,$ with $V_3 = -2c_1\ln S-c_2\ln E-c_3\ln I_m$, $V_4 = \left(S+E+I_m+I_s\right)^{m+1}$, $V_5 = -\ln S$, $V_6 = -\ln I_s,$ $c_i\in\mathbb{R}_+$ ($i = 1, 2, 3$); $M$ was a sufficiently large positive constant, while $m$ was a sufficiently small positive constant. Meanwhile, $M$ and $m$ satisfied the following relations:

where $B$, $F$, $e_1$, $e_2$ were defined in (4.6) and (4.12), respectively. Then, we obtained the following:

where $D_k = \left(\frac{1}{k}, k\right) \times \left(\frac{1}{k}, k\right) \times \left(\frac{1}{k}, k\right) \times \left(\frac{1}{k}, k\right).$ Obviously, $J\left(X\right)$ was a continuous function. We assumed that the minimum value of $J\left(X\right)$ was $\widetilde{J}$. We defined a non-negative $C^2$-function $Q\left(X\right) = J\left(X\right)-\widetilde{J}.$ Itô's formula acted on $V_3$, which gave the following:

By applying $a+b+c\geq3\sqrt[3]{abc}$ and $\sqrt[3]{a}+\sqrt[3]{b}\geq\sqrt[3]{a+b}$ for positive $a, b$ and $c$, expression (4.4) turned to the following:

We set

then,

Similarly,

with

Additionally, we derived the following:

From (4.5)–(4.7), we obtained the following:

We defined a bounded region as follows:

where $\varepsilon_1 > 0$ was sufficiently small and satisfied the following relations:

with

Obviously, $\mathbb{R}_+^4\backslash H = D_1\cup D_2\cup \cdots\cup D_{8}$, where

With (4.3) and (4.8), we discussed each case as follows:

Case 1. When $X\in D_1$, by (4.9), we obtained the following:

Case 2. When $X\in D_2$, due to (4.8), we derived the following:

Case 3. When $X\in D_3$, in view of (4.8), we obtained the following:

Case 4. When $X\in D_4$, according to (4.9), we derived the following:

Case 5. When $X\in D_5$, by (4.10), we obtained the following:

Case 6. When $X\in D_6$, due to (4.10), we derived the following:

Case 7. When $X\in D_7$, in view of (4.10), we obtained the following:

Case 8. When $X\in D_{8}$, according to (4.11), we derived the following:

Therefore, $\mathcal{L}Q\le-1$ when $X\in\mathbb{R}_+^4\backslash H$. Consequently, condition (ii) of Lemma 3.1 was satisfied. □

The sufficient conditions for the extinction were established with notation $\langle S(t)\rangle = \frac{1}{t}\int_{0}^{t}{S(s)}\text ds.$

Theorem 4.4. If the following conditions held,

then the solution of model (4.1) satisfies the following:

which implied that HIV/AIDS became extinct with an exponential rate.

Proof. Integrating the first equation of model (4.1), then, by Lemmas 4.1 and 4.2, we obtained the following:

which indicated

By Itô's formula, we defined $W\left(X\right) = E+I_m+I_s$ and obtained the following:

where

The integration on (4.14) provided the following:

Applying the strong law of large numbers ^[35], we obtained the following:

Recalling (4.13), (4.16), $R_0^e < 1$, and letting $t\rightarrow \infty$, we derived from (4.15) that

Thus, the number of infected individuals decreased to zero with an exponential rate in the long run. □

Remark 4.1. The following relationships held:

Moreover, $R_0 < R_0^e$ when $\min\{\gamma_e, \gamma_m\} \geq \frac{1}{3}\hat \sigma$. Precisely, model (4.1) had a unique stationary distribution when $R_0^s > 1$. Furthermore, $R_0 > R_0^s > 1$, which implied that the solution of model $(2.1)$ eventually approached the endemic equilibrium point $P^*$. However, $R_0 > 1$ revealed that HIV/AIDS became an endemic disease for model $(2.1)$, but not for model (4.1) , since $R_0^s > 1$ was not guaranteed.

5.   Parameter estimation and numerical simulations

5.1.   Features of surveillance data

In this section, we mainly analyzed the features of the surveillance data from the Fujian CDC in 2004–2011. Figure 1 presents the yearly incidence of HIV/AIDS infection cases with gender groups and age groups, and the maximum age, the minimum age and the average age for all HIV/AIDS infection cases in the Fujian Province are clearly presented within the left panel. Meanwhile, the yearly incidence of deaths with age groups are presented within the right panel of Figure 1, in which the cumulative incidence of deaths over eight years were 2.05/100,000 and 0.20/100,000 for all ages and for 60 years old and over, respectively. Precisely, the monthly incidence of infection cases and severe cases by genders were investigated in Figure 2.

5.2.   Estimation for parameters and optimal simulations

According to the report in ^[40], the total population scale and the life expectancy of the Fujian Province in 2004 were 35,290,000 and 73.834, respectively, which then gave the following initial values of this study:

As of January 2004, nine HIV/AIDS infection cases were recorded by the Fujian CDC, in which seven mild cases (72.7%) and two severe cases (27.3%) were reported with no deaths. According to the research results in ^[41], the basic reproduction number of HIV/AIDS ranged from 1.3 to 6.0. Subsequently, we computed the total HIV/AIDS infection cases as follows: $7\times 6.0 + 2\times 1.3 = 44.6\approx 45.$ By proportion, we derived $I_m(0) = 33$ and $I_s(0) = 12$ because we took potential HIV/AIDS infection cases before 2004 into account in this study. Furthermore, we assumed that the number of HIV/AIDS infection cases with no clinical symptoms was 1,000, that is, $E(0) = 1,000$. By $S(0) = N(0)-E(0)-I_m(0)-I_s(0)$, we computed that $S(0) = 35,288,955$. Therefore, all initial values were given herewith.

With modified HIV policies from the Chinese government in February 2006 in ^[42], the targeted population with HIV test was extended, which produced an increase in the number of HIV/AIDS infection cases from 2006. Therefore, the values of $\beta_e$ were usually regarded as distinct before and after 2006. Utilizing the least squares method, $\beta_e$ was estimated by two parts: $\beta_e^1 = \frac{0.258}{35,290,000}$ for 2004–2005, $\beta_e^2 = \frac{0.417}{35,290,000}$ for 2006–2011. Moreover, $\beta_m$ was estimated by $\frac{0.120}{35,290,000}$. By comparing the population scales of Mexico and the Fujian Province, we took $\gamma_e = 0.1411$ in ^[43]. The surveillance data from the Fujian CDC indicated that the probability that an HIV/AIDS infection case developed into a severe case was $0.2731$. Meanwhile, the average time from $I_m$ to $I_s$ was $2.0162$ years, which gave $\gamma_m = \frac{0.2731}{2.0162}$. The average time for a severe case to death was $1.1419$ years, which indicated that $\delta = \frac{1}{1.1419}$. We collected all parameters in Table 1.

By the surveillance data from the Fujian CDC, the gender and age distributions of HIV/AIDS infection cases in 2004–2011 were recorded, in which the proportions for male infection cases (69.9%) and female infection cases (30.1%) were investigated during the period 2004–2011, and the percentages for key infection cases (85.9%, 20–49 years old) and ordinary infection cases (14.1%, 1–19 years old as well as 50 years old and over) were discussed by clustering the data. The curves for cumulative incidences of infection cases by genders are plotted on left panel of Figure 3, in which the simulations for total incidence (in yellow), male incidence (in blue) and female incidence (in pink) were performed. The analysis for the cumulative incidences of infection cases by age are shown on right panel of Figure 3. By the same argument, similar investigations on the cumulative incidence of severe cases, the proportions of male and female cases, and the proportions of key and ordinary cases are presented in Figure 4.

Furthermore, predictions of HIV/AIDS incidence for the period 2012–2014 were made. We chose the values of $\beta_e$ from $\frac{0.377}{35,290,000}$ to $\frac{0.417}{35,290,000}$, and kept other parameters the same as Table 1. Then, the results of the predictions were derived, which ranged from 22.97/100,000 to 23.84/100,000 at the end of 2014 as presented on left panel of Figure 5. The minimum value of $\beta_e$ was estimated by $\frac{2\beta_e^1+6\beta_e^2}{8}$, and the maximum value was same with $\beta_e^2$. By the proportion of male infection cases, the possible incidence rate ranged from 16.6/100,000 and 16.66/100,000, and by the percentage of key infection cases, the possible incidence rate ranged from 19.73/100,000 to 20.48/100,000 at the end of 2014. Especially, by the same discussion, the predictions of incidence for severe cases during the period 2012–2014 were operated, which produced the range from 7.12/100,000 to 7.23/100,000 at the end of 2014, and were presented on right panel of Figure 5.

5.3.   Comparisons for dynamical properties

The recent contributions in ^[44,45] established approaches regarding the positivity preserving truncated Euler-Maruyama (PPTEM) method, which supported the numerical simulations on the extinction and persistence of stochastic models. Moreover, the PPTEM method was governed on the extinction of stochastic models to avoid the negative values. We used the Milstein's higher order (MHO) method ^[46] to establish the discretization equations of model (4.1), due to its efficiency and performance of the persistence of stochastic models.

As given by Table 1, the basic reproduction number $R_0$ was 2.4032 for the period between 2004 and 2005 and 3.3990 for the period between 2006 and 2011. Our investigation showed that the solution of model (2.1) eventually approached the globally asymptotically stable $P^*$, but was repelled by the unstable $P_0$ as claimed in Theorem 3, where $P^*$ and $P_0$ were calculated as $P^* = (2.91 \times 10^4, 6.21 \times 10^3, 5.88 \times 10^3, 8.59 \times 10^2)^\text{ T},$ and $P_0 = (10^5, 0, 0, 0)^\text{ T}.$ Additionally, our investigation revealed that model (4.1) admitted the indicator $R_0^s \approx 3.3849$ when taking $\sigma_i = 0.01 (i = 1, 2, 3, 4)$. By Theorem 4.3, we concluded that HIV/AIDS was stochastically persistent for a long run as demonstrated on left panel of Figures 6–8. When we repeated 3,000 runs, the corresponding curves for the histogram were respectively carried out as shown on right panel of Figures 6-8. Specifically, we noticed that $R_0 > R_0^s > 1$ held in Remark 4, which indicated that the stochastic persistence to model (4.1) implied the persistence to model (2.1) for HIV/AIDS.

Furthermore, we performed an operation for models (2.1) and (4.1) in $400$ years. The research results showed that the solution for model (4.1) and the solution for model (2.1) matched well when we took $\sigma_i = 0.01(i = 1, 2, 3, 4)$. The deviations are shown in Figure 9. Then, we defined the proportion of deviation for the solutions of models (2.1) and (4.1) by the following form:

where the proportion of deviations varied with the time $t$. The research results showed that when we took $\sigma_i = 0.01(i = 1, 2, 3, 4)$, the maximum absolute values of the proportion of deviations were $5.35\%$ for $E$, $5.93\%$ for $I_m$ and $6.15\%$ for $I_s$.

Furthermore, the changes of the solutions and the probability density functions to model (4.1) were explored when the intensities of white noises $\sigma_i = 0.01, 0.02, 0.03, 0.04, 0.05$ $(i = 1, 2, 3, 4)$ were chosen and other parameters were kept the same as Table $1$. Then, the corresponding numerical simulations of the solutions were derived on left panels of Figures 10–12. Meanwhile, the changes of the probability density functions were demonstrated on right panels of Figures 10–12 when the intensities of white noises increased.

The research results showed that the tendencies of the numerical simulations for next five decades were very similar when we compared Figure 6.1 for Indonesia in ^[25] and Figures 7–8 for the Fujian Province in this study. The incidence (1/100,000) of HIV/AIDS infection cases for the Fujian Province reached a peak with about $2.86 \times 10^4$, then rapidly declined and gradually became stable around at $1.61 \times 10^4$, which revealed that the potential risk of medical burden raised near the infection peak. We suggested that the local government adopted more tough strategies before the medical burden appeared, and medical service centers were asked to prepare the future medical resources as early as possible.

6.   Conclusions

By surveillance data from 2004 to 2011 of the Fujian CDC, we proposed symptom-dependent HIV/AIDS models since the mild and severe symptom cases were recorded well, and further investigated the survival dynamics such as persistence, extinction and stationary distribution. The research results showed that the disease-free equilibrium $P_0$ was locally and globally asymptotically stable for model (2.1) when $R_0 < 1$; additionally, $P_0$ was unstable, whereas the endemic equilibrium point $P^*$ was locally and globally asymptotically stable when $R_0 > 1$. When it came to the environmental fluctuations, the interaction mechanism among compartments of model (2.1) remained the same. Based on the existence and uniqueness of the global positive solution, the ergodic stationary distribution of stochastic model (4.1) was valid when $R^s_0 > 1$, which indicated that HIV/AIDS prevailed for the long run. Moreover, we investigated the sufficient condition of extinction for HIV/AIDS for a short period of time.

The surveillance data of HIV/AIDS from the Fujian CDC gave that $R_0 \approx 3.3990 > 1$ to model (2.1), which indicated that HIV/AIDS became an endemic disease for a long time in the Fujian Province, and that HIV/AIDS displayed an exponential upward trend in Figure 5. Meanwhile, model (4.1) displayed small fluctuations (say $\sigma_i = 0.01$) provided that $R_0^s \approx 3.3849 > 1$, which revealed that HIV/AIDS cases of the Fujian Province persisted almost surely, as shown in Figures 6-8. However, when $\sigma_i = 0.10$ were taken in Theorem 4 such that $R_0^s \approx 2.3876 > 1$, the sample paths did not present the stationary properties very well. To perform the 90–90–90% plan of WHO in ^[47], the scale to the HIV infection mainly relied on the reductions of $\beta_e$ (the effective contact coefficient between $S$ and $E$) and $\beta_m$ (the effective contact coefficient between $S$ and $I_m$). For the susceptible population, the reduction of HIV infection might come from the following: individual-oriented protection such as male or female condom use; medical-system-oriented protection such as HIV tests and voluntary medical male circumcision in ^[48]; and social-oriented protection such as harm reduction services for people who inject and use drugs. Alternatively, some medicines and medical devices such as antiretroviral drugs (ARVs), dapivirine vaginal rings and injectable long acting cabotegravir were extensively available to prevent HIV, which reduced the value of $\gamma_m$, and the numbers of severe cases and deaths were therefore declined.

Figure 5 indicated that the HIV/AIDS incidence rate ranged from 22.97/100,000 to 23.84/100,000 in the Fujian Province at the end of 2014; in comparison, the research results in ^[25] showed that the HIV/AIDS incidence rate was 36.36/100,000. In 2020, the incidence rate of HIV/AIDS infection cases for Mainland of China admitted 76.47/100,000 in ^[25], while the prediction by model (2.1) showed that the lower and upper incidence rates for the Fujian Province were 120.16/100,000 and 149.04/100,000, respectively, at the end of 2020. In comparison with other previous contributions, our models consisted of four equations with seven parameters to handle the surveillance data from the Fujian CDC, and the optimal fittings were performed very well. Our model (2.1) modified model (2.1) in ^[25] by combining the equations for $I$ and $C$ as an equation of $I_m$, and keeping the equation of $A$ as an equation of $I_s$.

Use of AI tools declaration

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

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgment

Authorship contribution statement: Wenshuang Li: Modelling, Data curation, Software, Visualization, Writing-original draft. Shaojian Cai: Data curation, Validation, Writing-review. Xuanpei Zhai: Writing-original draft. Jianming Ou: Writing-review. Kuicheng Zheng: Writing-review. Fengying Wei: Methodology, Validation, Writing-review & editing. Xuerong Mao: Supervision, Methodology, Writing-review.

The research of F.Wei was supported in part by Natural Science Foundation of Fujian Province of China (2021J01621); the research of X. Mao was supported by the Royal Society, UK (WM160014, Royal Society Wolfson Research Merit Award), the Royal Society of Edinburgh (RSE1832), Engineering and Physical Sciences Research Council, UK (EP/W522521/1) for their financial support.