A two-group epidemic model with heterogeneity in cognitive effects

1.   Introduction

In recent years, infectious diseases such as Dengue Fever (1979), severe acute respiratory syndromes (SARS) (2003), influenza a (H1N1) (2009), and Corona virus disease 2019 (COVID-19) (2019) have spread globally. The impact of these infectious diseases extends beyond the health sector, severely threatening the lives and health of people worldwide ^[1,2,3,4,5]. In the age of information, the spread of infectious diseases is often accompanied by the rapid dissemination of information. On the one hand, information related to diseases can raise public awareness, providing details about transmission routes, infectiousness, and possible preventive measures, thereby enabling individuals to take effective protective actions. On the other hand, disease-related information can also cause panic, leading to negative behaviors. As a result, the spread of disease-related information can induce changes in individual behavior, which, in turn, significantly impacts the transmission of the disease ^[6,7,8].

With the continuous rise of social media, news channels, scientific research, and other information sources, public understanding and responses to infectious diseases are undergoing profound changes. Media information not only serves as a carrier of knowledge, but also plays a pivotal role in disease prevention and control. Many ordinary differential equation models were used to analyze the impact of individual behavioral changes, such as wearing masks, maintaining social distancing, etc., in response to the spread of disease-related information. Generally, studies examining the dissemination of disease-related information can be classified into two main approaches: one assumes that behavioral changes lead to a reduction in infection or contact rates. This approach typically models infection rates as a function of the number of infected individuals ^[9,10,11,12] or the number of news items ^{[13,14,15,16]}, treating news items as a separate compartment. For example, Luo and Liu ^[17] studied an susceptible-vaccinated susceptible-exposed-infected-quarantined-recovered (SSvEIQR) model with nonlinear contact rate, isolation rate and vaccination rate driven by media coverage. The sufficient conditions to prove the global stability of the endemic equilibrium were obtained by applying the geometric method into the four-dimensional system. They obtained some measures to control the spread of the disease, such as reducing contact, strengthening isolation, and vaccination. The second approach divides the population into two categories: those who are sensitive to disease-related information and those who are not. Individuals sensitive to such information may subconsciously take protective measures, thus reducing exposure rates. This method introduces a heterogeneous modeling solution, which inevitably increases the system's dimensionality and complicates theoretical analysis ^[18,19,20]. For instance, Li and Xiao ^[21] proposed a two-group model that simulates both disease and information propagation, incorporating saturated recovery rates. By applying fast-slow theory and analyzing the system dynamics directly, they identified the existence and stability of potential equilibrium points and the occurrence of backward bifurcation. Consequently, considering heterogeneous models becomes crucial, as variations in information sensitivity present significant challenges in controlling and eliminating infectious diseases.

In the early stages of emerging infectious diseases, the most significant challenge is the limitation of medical resources, a threat faced by nearly all countries. For instance, during the initial outbreak of COVID-19 in mainland China in early 2020, the shortage of detection kits, hospital beds, and ventilators emerged as the most pressing issue. The need to ensure an adequate supply of medical resources to effectively respond to outbreaks has become increasingly emphasized. In this context, Cui et al. ^[22,23] developed an epidemiological model incorporating saturated recovery in infected individuals, demonstrating that saturated recovery can lead to bistability and periodicity. To analyze the impact of medical conditions, Wang ^[24] proposed an Ebola epidemic model that incorporates limited medical resources, immunity loss, and the tracking and quarantining of susceptible individuals. The study explores bifurcation phenomena and conducts sensitivity analysis, revealing a strong correlation between the control reproduction number and the incineration or burial rate of dead bodies. Asamoah et al. ^[25] introduced a mathematical model to describe the nonlinear recovery rate in bacterial meningitis, offering a framework for controlling disease transmission in resource-limited settings, and identifying both forward and backward bifurcations. Li and Xiao ^[26] proposed an epidemic $SEIM$ model incorporating saturated media growth and a saturated recovery rate. The theoretical results suggest that only nonlinear recovery models may exhibit backward bifurcation under specific conditions. Through numerical simulations, their model revealed a rich array of dynamical behaviors, including forward and backward bifurcations, Hopf bifurcations, saddle-node bifurcations, homoclinic bifurcations, and unstable limit cycles. Thus, the limitations of medical resources give rise to complex dynamics, posing significant challenges in the elimination and control of infectious diseases.

As discussed in ^[21,26,27], during the outbreak of emerging infectious diseases, medical resources are often in short supply, and leveraging media information to control disease transmission proves to be a highly effective strategy. Consequently, it is essential to incorporate the limitations of medical resources into epidemiological modeling. While existing research on mathematical models incorporating media information is relatively well-developed, few studies have considered the heterogeneity of different population groups or examined how media-influenced infection rates and limited medical resources interact to impact disease dynamics. Furthermore, the challenge of formulating more effective and targeted control strategies under conditions of limited media resources remains an area of concern. In light of these issues, this paper aims to explore the heterogeneity in population responses to media and the constraints posed by limited medical resources. Specifically, we focus on how these two factors influence disease transmission and control strategies within a two-group model, taking into account the limitations of media resources.

The structure of this paper is organized as follows. Section 2 presents the model formulation for heterogeneous infectious diseases. In Section 3, we analyze the well-posedness of the model, calculate the basic reproduction number, and investigate the local stability of the disease-free equilibrium. Additionally, we establish the conditions for disease extinction and persistence in the absence of medical resource constraints. In Section 4, we formulate a control problem that considers media information allocation and solve the corresponding optimal control problem with control variables. Section 5 is dedicated to numerical simulations, where we explore the sensitivity of the model to parameter changes, the bifurcation structure, and the effectiveness of two types of control strategies. Finally, Section 6 provides a summary of the findings and an outlook for future research.

2.   Model formulation

We use a classical SIR model to describe the dynamics of disease transmission, where $S$, $I$, and $R$ represent the susceptible, infectious, and removed (recovered or deceased) individuals, respectively. To account for individual responses to the information they receive, we divide the total population into two groups: information-sensitive groups ($S_1$, $I_1$, $R_1$) and information-insensitive groups ($S_2$, $I_2$, $R_2$) ^[21]. We assume that the transmission rate between group $j$ and group $i$ is given by $\beta_{ij}$, where $i, j = 1, 2$. Due to the limitations of medical resources, infectious individuals enter the removed compartment with a recovery rate that follows a saturation function, $\frac{\gamma_i}{1 + h_i I_i}$, where $i = 1, 2$. Here, $\Lambda_i$ (for $i = 1, 2$) represents the immigration rates, and $d$ is the natural mortality rate for all individuals.

To simulate the impact of media reports and information on disease transmission, we consider the number of media reports related to the epidemic as an independent variable, denoted as $M(t)$. To capture the effect of media reports in reducing the effective contact between susceptible and infected individuals, we introduce an exponential decay factor, as proposed in ^[26], represented by $e^{-\alpha_1 M}$ and $e^{-\alpha_2 M}$, which describes the declining coefficient of incidence for the susceptible and infected groups, respectively. Parameters $\mu_1$ and $\mu_2$ represent the response rates of infected individuals who are sensitive and insensitive to information, respectively. Finally, $\tau$ denotes the natural disappearance rate of media reports. Based on this framework, we consider the following model of differential equations:

The model flow diagram is shown in Figure 1, and all parameters are nonnegative constants, and their descriptions, values, and sources are provided in Table 1. Since the equations for $R_1$ and $R_2$ can be decoupled from the other equations in model (2.1), we exclude these equations from further consideration in this work. Therefore, we focus on discussing the following simplified model,

3.   Model analysis

3.1.   Well-posedness

We first prove the nonnegativity and boundedness of the solutions of model (2.2) to show the the biological feasibility of it.

Theorem 3.1. Model (2.2) has a bounded solution for all time $t \ge 0$ with initial condition lying in domain

Moreover, the compact set $\Omega$ is positively invariant with respect to model (2.2).

Proof. Since the vector field described by the righthand sides of model (2.2) is Lipschitz continuous in $\Omega$, a unique solution exists when $t \ge 0$. Note that $S_1 = 0 \Rightarrow S_1' \ge 0$, $I_1 = 0 \Rightarrow I_1' \ge 0$, $S_2 = 0 \Rightarrow S_2' \ge 0$, $I_2 = 0 \Rightarrow I_2' \ge 0$, and $M = 0 \Rightarrow M' \ge 0$. Therefore, all solutions of model (2.2) starting from $\Omega$ remain in $\Omega$, and model (2.2) is mathematically and epidemiological well-posed in $\Omega$.

The total population is represented by $N = S + I + R$. Depending on the difference of groups, the total population $N$ is divided into two parts, $N_1$ and $N_2$, to represent the total population number of two groups. To demonstrate the boundedness of the solutions of model (2.1), we add all equations of model (2.1):

The above inequality implies that $N$ is bounded above as well as below. Now we prove the boundedness of $M$ through the equation of $M$ in (2.1). By some calculations as,

Thus $0 \le M \le \frac{{{\mu _1}{\Lambda _1} + {\mu _2}{\Lambda _2}}}{{\tau d}}$, which indicates that the feasible region for the model (2.2) is

□

3.2.   Disease-free equilibrium and basic reproduction number

We now demonstrate that model (2.2) has a disease-free equilibrium (DFE) given by

Additionally, we determine the basic reproduction number ($\mathcal{R}_0$), which serves as a critical threshold in epidemiological models. $\mathcal{R}_0$ represents the average number of secondary infections generated by an infectious individual in a population of susceptible individuals during the infectious period. To calculate $\mathcal{R}_0$, we use the next-generation matrix approach ^[28] and rewrite the equations for the infectious compartments as follows:

and

Then we calculate the Jacobian matrices $F$ and $V$ at the $E_0$, which are given as

The basic reproduction number is derived by computing the spectral radius of the next-generation matrix, $F V^{-1}$,

where

Next, we examine the dynamical behavior of model (2.2) around the DFE $E_0$, and identify the basic reproduction number $\mathcal{R}_0$ as a crucial threshold parameter that governs the disease's extinction or persistence. To begin, we analyze the local stability of the $E_0$, by evaluating the sign of the eigenvalues of the Jacobian matrix at the $E_0$.

Theorem 3.2. The DFE $E_0$ of model (2.2) is locally asymptotically stable (LAS) when $\mathcal{R}_0 < 1$, and it becomes unstable when $\mathcal{R}_0 > 1$.

Proof. The Jacobian matrix at $E_0$ is given by:

It is clear that the eigenvalues $-d$, $-d$, and $-\tau$ are three negative eigenvalues of $J(E_0)$. Therefore, the local stability of the DFE, $E_0$, depends on the remaining eigenvalue of $J(E_0)$. The sign of the remaining eigenvalues is determined by the rest of the Jacobian matrix, denoted as $J_1(E_0)$. The eigenvalues of $J_1(E_0)$ are obtained by solving the characteristic equation:

Thus, the characteristic equation is:

In the light of Routh-Hurwitz criterion, we need to prove

and

When $\mathcal{R}_0$ is less than 1, we can get ${k_1} + {k_2} + \sqrt {{{({k_1} - {k_2})}^2} + 4{k_3}{k_4}} < 2$. Upon further simplification, this leads to the condition

After substitution, we can get

which yields

Thus, we have $\lambda_1 \lambda_2 > 0$. Next, we demonstrate that $\lambda_1 + \lambda_2 < 0$. Since $\mathcal{R}_0 > 0$, we deduce that $k_1 + k_2 + \sqrt{(k_1 - k_2)^2 + 4 k_3 k_4} > 0,$ which implies that $k_3 k_4 > k_1 k_2$. Substituting into the specific expressions, we obtain $\beta_{12} \beta_{21} > \beta_{11} \beta_{22}.$ Furthermore, from the inequality

we substitute $\beta_{12} \beta_{21} > \beta_{11} \beta_{22}$ into the above inequality to obtain the result:

After simplification, we get $\frac{{{\beta _{11}}{\Lambda _1}}}{{d\left({d + {\gamma _1}} \right)}} + \frac{{{\beta _{22}}{\Lambda _2}}}{{d\left({d + {\gamma _2}} \right)}} < 1$, which implies $\frac{{{\beta _{11}}{\Lambda _1}}}{{d\left({d + {\gamma _1}} \right)}} < 1$ and $\frac{{{\beta _{22}}{\Lambda _2}}}{{d\left({d + {\gamma _2}} \right)}} < 1$, and this can be further written as $\beta_{11} \Lambda_1 < d(d + \gamma_1)$ and $\beta_{22} \Lambda_2 < d(d + \gamma_2).$ Therefore, we have

Through the calculation, when $\mathcal{R}_0$ is less than 1, we can get $\text{tr}\left(J(E_0) \right) < 0$ and $\text{det}\left({J({E_0})} \right) > 0$. Therefore, according to the Routh-Hurwitz criterion, we complete the proof. □

3.3.   Dynamics of the model (2.2) without medical resource limitation

Note that model (2.2) has high nonlinearity, which makes the analysis of the model's dynamics challenging. Therefore, we consider a special case:

This special case only includes the factor of contact rate influenced by information, i.e., $h_1 = 0$ and $h_2 = 0$. Then, we proceed to theoretically analyze the conditions for the extinction or persistence of the disease in the special case.

Theorem 3.3. Assume model (3.1) satisfies $\mathcal{R}_0 < 1$, then $E_0$ is globally asymptotically stable in $\Omega$ given there exist two positive constant $\rho _1$ and $\rho _2$ with

and

Proof. In model (3.1), we can continue to use the conditions of the Routh-Hurwitz criterion in Section 3.2, because when $h_1$ = 0 and $h_2$ = 0, the conditions are interlinked. We use the method of constructing the Lyapunov function to prove the disease extinction of model (3.1). The Lyapunov function is given by:

It is straightforward that $L$ is a positive definite function. Along the trajectories of model (3.1), we have

Note that the first and second terms are clearly nonpositive, and the sign of the derivative of the Lyapunov function is primarily determined by the latter two terms. Specifically, we focus on the expression inside the parentheses of the third term,

Therefore, if $\frac{{{\beta _{11}}{\Lambda _1}}}{{d\left({{\gamma _1} + d} \right)}} \le \frac{1}{2}$ and $\frac{{{\beta _{21}}{\Lambda _2}}}{{d\left({{\gamma _1} + d} \right)}} \le \frac{1}{2}$ are satisfied, it can be demonstrated that the third term is nonpositive. To achieve this, we assume that there exists $0 < {\rho _1} < 1,$ such that $0 < \frac{{{\beta _{11}}{\Lambda _1}}}{{d\left({{\gamma _1} + d} \right)}} < {\rho _1}.$ Due to $\frac{{{\beta _{21}}{\Lambda _2}}}{{d\left({{\gamma _1} + d} \right)}} = \frac{{{\beta _{11}}{\Lambda _1}}}{{d\left({{\gamma _1} + d} \right)}} \cdot \frac{{{\beta _{21}}{\Lambda _2}}}{{{\beta _{11}}{\Lambda _1}}},$ we can rewrite the expression as $\frac{{{\beta _{21}}{\Lambda _2}}}{{d\left({{\gamma _1} + d} \right)}} = \frac{{{\beta _{11}}{\Lambda _1}}}{{d\left({{\gamma _1} + d} \right)}} \cdot \frac{{{\beta _{21}}{\Lambda _2}}}{{{\beta _{11}}{\Lambda _1}}} < {\rho _1} \cdot \frac{{{\beta _{21}}{\Lambda _2}}}{{{\beta _{11}}{\Lambda _1}}}.$ To ensure that $\frac{{{\beta _{21}}{\Lambda _2}}}{{d\left({{\gamma _1} + d} \right)}} \le \frac{1}{2}$, it is sufficient to have ${\rho _1} \cdot \frac{{{\beta _{21}}{\Lambda _2}}}{{{\beta _{11}}{\Lambda _1}}} < \frac{1}{2},$ which implies ${\rho _1} < \frac{{{\beta _{11}}{\Lambda _1}}}{{2{\beta _{21}}{\Lambda _2}}}.$ Thus, it is enough to require that $\max \left\{ {\frac{{{\beta _{11}}{\Lambda _1}}}{{d({\gamma _1} + d)}}, \frac{{{\beta _{21}}{\Lambda _2}}}{{d({\gamma _1} + d)}}} \right\} < {\rho _1} < 1$ in order to guarantee that the third term is nonpositive. The treatment of the last term follows a similar approach. We can conclude that there exists $0 < {\rho _2} < 1,$ and it is sufficient to require that ${\rho _2} < \frac{{{\beta _{22}}{\Lambda _2}}}{{2{\beta _{12}}{\Lambda _1}}}$ and $\max \left\{ {\frac{{{\beta _{12}}{\Lambda _1}}}{{d({\gamma _2} + d)}}, \frac{{{\beta _{22}}{\Lambda _2}}}{{d({\gamma _2} + d)}}} \right\} < {\rho _2} < 1.$ These conditions ensure that the last term is also nonpositive.

Hence, together with conditions, we have $\frac{{dL}}{{dt}} \le 0$. Then, $\frac{{dL}}{{dt}} = 0$, if ${S_1} = S_1^0, {I_1} = I_1^0, {S_2} = S_2^0, {I_2} = I_2^0.$ Therefore, the largest invariant set contained in $\left\{ {\left({{S_1}, {I_1}, {S_2}, {I_2}} \right)\left| {\frac{{dL}}{{dt}} = 0} \right.} \right\}$ is $\left({S_1^0, I_1^0, S_2^0, I_2^0} \right).$ According to the Lasalle's invariance principle ^[29], the DFE is globally asymptotically stable. □

When $\mathcal{R}_0 > 1$, $E_0$ becomes unstable. Therefore, it is reasonable to assume that in this case, the infected population ${I_1}\left(t \right)$ and ${I_2}\left(t \right)$ will remain continuous. In fact, the following theorem holds.

Theorem 3.4. Model (3.1) is uniformly persistent when $\mathcal{R}_0 > 1$ if, and only if, there exists a nonnegative constant $\eta > 0$, such that all solutions of model (3.1) satisfy:

Proof. We will apply Theorem 3.4 in ^[30] to prove the uniform persistence of model (3.1). Set

which is relatively closed in $X$.

Let $\Phi \left(t \right):X \to X$ be the solution flow associated with model (3.1), that is, $\Phi \left(t \right)\left({{A_0}} \right) = A\left(t \right)$. According to the positivity, it is easy to know when ${S_1}\left(0 \right) \ge 0, {I_1}\left(0 \right) \ge 0, {S_2}\left(0 \right) \ge 0, {I_2}\left(0 \right) \ge 0, M\left(0 \right) \ge 0$, $\Phi \left(t \right)$ is positively invariant about $X$, so for $\forall t > 0$, then ${S_1}\left(t \right) > 0, {I_1}\left(t \right) > 0, {S_2}\left(t \right) > 0, {I_2}\left(t \right) > 0, M\left(t \right) > 0$. Now we prove that ${X_0}$ is positive invariant for $\Phi \left(t \right)$. According to the second and fourth equations of model (3.1), we get that

If ${I_1}(0) \ge 0, {I_2}(0) \ge 0$, according to the above two equations, we can get

Therefore, ${X_0}$ is positive invariant for $\Phi \left(t \right)$. Moreover, model (3.1) has the ultimate upper boundary, thus we obtain $X$ is point dissipative for $\Phi \left(t \right)$. We set ${M_\partial } = \left\{ {{A_0} \in \partial {X_0}, \forall t \ge 0} \right\}.$ We claim that ${M_\partial } = \left\{ {\left({{S_1}, 0, {S_2}, 0, 0} \right)} \right\}.$

Hypothesize that $A\left(t \right) \in {M_\partial }, \forall t \ge 0$, and it suffices to show that ${I_1}(t) = {I_2}(t) = 0, \forall t \ge 0$. If it is not true, then there exists ${t_0} \ge 0$ such that either (a) ${I_1}({t_0}) > 0, {I_2}({t_0}) = 0$; or (b) ${I_1}({t_0}) = 0, {I_2}({t_0}) > 0$.

For case (a), from the fourth equation of model (3.1), we have

Hence, there is an ${\varepsilon _0} > 0$ such that ${I_2}(t) > 0$, $\forall t \in \left({{t_0}, {t_0} + {\varepsilon _0}} \right)$. On the other hand, from ${I_1}\left({{t_0}} \right) > 0$, there exists an $\varepsilon _1$ $\left({0 < {\varepsilon _1} < {\varepsilon _0}} \right)$ such that ${I_1}\left(t \right) > 0$, $\forall t \in \left({{t_0}, {t_0} + {\varepsilon _1}} \right)$. Thus, we have ${I_1}\left(t \right) > 0, {I_2}\left(t \right) > 0, \forall t \in \left({{t_0}, {t_0} + {\varepsilon _1}} \right)$, which contradicts the assumption that $\left\{ {{S_1}\left(t \right), {I_1}\left(t \right), {S_2}\left(t \right), {I_2}\left(t \right), M\left(t \right) \in {M_\partial }, \forall t \ge 0} \right\}$. Similarly, we can obtain a contradiction for case (b). This proves the claim.

Let $\mathcal{A} = { \cap _{x \in {A_b}}}\omega \left(x \right)$, where ${{A_b}}$ is the global attractor of model (3.1) restricted to $\partial {X_0}$. We show that $\mathcal{A} = {E_0}$. In fact, from $\mathcal{A} \subseteq {M_\partial }$ and the first and third equations of model (3.1), we have

Thus, $E_0$ is the isolated invariant set in $X$.

Next we show that ${W^s}\left({{E_0}} \right) \cap {X_0} = \emptyset$. If it is not true, then there exists a solution $\left({{S_1}(t), {I_1}(t), {S_2}(t), {I_2}(t), M(t)} \right) \in {X_0}$ such that

Therefore, for any sufficiently small constant $\xi > 0$, there exists a positive constant ${T_1} = {T_1}\left(\xi \right)$ such that we have the following inequalities for all $t \ge {T_1}$:

Consider the following auxiliary system with the above constant $\xi$:

Considering the following comparing system, where ${{\tilde I}_1}\left(0 \right) = {I_1}\left(0 \right), {{\tilde I}_2}\left(0 \right) = {I_2}\left(0 \right):$

Then, the nonnegative matrix $\tilde F$ and the non-singular M-matrix $\tilde M$, represented as new-infection and transition matrices severally, for the model (3.1), are provided by

Denote

According to Lemma 2.1 in ^[31], as long as $\mathcal{R}_0 > 1$, when $\xi$ and $s\left({\tilde J\left(0 \right)} \right) > 0$, $s\left({\tilde J\left(\xi \right)} \right)$ is continuous. We can select such sufficiently small $\xi > 0$ that we obtain $s\left({\tilde J\left(\xi \right)} \right) > 0$, and this can explain that the positive solutions of the lower comparing system increase exponentially. By contrasting to the standard comparing system, as $t \to \infty$, the solutions ${I_1}\left(t \right)$ and ${I_2}\left(t \right)$ of model (3.1) are infinite, which are opposite to the truth that the solutions of the model (3.1) are bounded. Thus, ${W^s}\left({{E_0}} \right) \cap {X_0} = \emptyset$. Clearly, every orbit in ${M_\partial }$ converges to ${E_0}$. Thus, by Theorem 3 in ^[32], we have

By Theorem 4.3 and Remark 4.3 in ^[33], we conclude that model (3.1) is uniformly persistent with respect to $\left({{X_0}, \partial {X_0}} \right)$. From Theorem 2.4 in ^[34], model (3.1) has at least one equilibrium $\left({S_1^*, I_1^*, S_2^*, I_2^*, {M^*}} \right) \in {X_0}$, with $I_1^* \gg 0$ and $I_2^* \gg 0$. Accordingly, there is at least one an endemic equilibrium in model (3.1). The proof of the theorem is complete. □

4.   Optimal control under limited media resources

In this section, we use optimal control techniques to study the model (2.1). In order to minimize the cost of implementing control strategies, it is necessary to find time-dependent control strategies. Most of the control strategies used in daily life are considering continuous control. In fact, this problem is a typical optimal control problem. If public health institutions want to eradicate a disease, maintaining a high level of control is crucial. However, this often incurs significant economic costs, so it is necessary to develop a time-dependent control strategy. The control strategy we adopt focuses on the allocation of media resources. Specifically, the total media resource allocation for two groups of infected individuals sums to 1. Let the media resource allocation rate for the sensitive infected individuals be denoted as $u(t)$, while the allocation rate for the insensitive infected individuals is $1 - u(t)$. Our primary objective is to minimize both the total number of infections and the costs associated with adjusting media coverage intensity.

Another key issue is the limitation of media resources when using media information to control infectious diseases. In particular, during emergencies, media resources are often constrained, as platforms such as television, radio, and social media have limitations in terms of coverage area, transmission frequency, and target audience. Therefore, communicators must prioritize the delivery of essential information and avoid excessive or redundant details. In the optimal control problem we consider, we assume a limited media resource, denoted as $V$. Optimal control problem with minimizing objective function

subject to

where the initial condition is the same as that of model (2.1). The coefficients $A$, $B$ and $C/2$ are positive. $A$, $B$ represent the weight of the two groups of infected people respectively. Here we assume that $C/2$ is the weight associated with control ${u(t)}$. Note that ${u(t)}$ is a Lebesgue measurable function on a finite interval $\left[{0, {t_f}} \right]$, where $0 < u(t) < 1$. To start, we prove the existence of an optimal control for the model (4).

Theorem 4.1. There exists an optimal control ${u^*}(t)$ such that

subject to the control model (4.1).

Proof. By the result in ^[35], we prove the existence of an optimal control. Note that the control and the state variable are nonnegative values. In this minimizing problem, the necessary convexity of the objective functional in ${u(t)}$ is satisfied. Meanwhile, ${u(t)}$ belongs to the control set $U$, where

The boundedness of the optimal system ensures the compactness required for the existence of an optimal control. In addition, the integrand in the objective function is convex on the control set $U$. Also, we can see that there exist a constant $\rho > 1$ and positive numbers ${\omega _1}, {\omega _2}$ such that

because the state variables are bounded, which completes the existence of an optimal control. □

In order to find the optimal solution, we use Pontryagin's maximum principle as follows ^[36]. Define the Hamiltonian $H$ for the control problem:

Then, the adjoint equations are given:

By the optimality conditions, we have

Note that the boundness is placed on the control variable ${u(t)}$, then the optimality condition is changed to

which is the optimal control.

5.   Numerical simulations

5.1.   Influence of media and limited medical resources on the number of infected individuals

First, we study the impact of the reaction of susceptible individuals to the media, $\alpha_1$ and $\alpha_2$, on the number of sensitive infected individuals, insensitive infected individuals, and total infected individuals at steady state, which are given in Figure 2. Figure 2(a) shows that $I_1$ decreases with the increase of $\alpha_1$, however, the impact of $\alpha_2$ remains negligible. Figure 2(b) illustrates that a higher value of $\alpha_2$ leads to a more significant reduction in the number of insensitive infected individuals. Figure 2(c) indicates that increases in both $\alpha_1$ and $\alpha_2$ contribute to a decrease in the total number of infected individuals. Figure 3 shows the impact of the response rate of infected individuals, $\mu_1$ and $\mu_2$, on the number of sensitive infected individuals, insensitive infected individuals, and total infected individuals at steady state. As shown in Figure 3(a), both an increase in $\mu_1$ and $\mu_2$ lead to a decrease in $I_1$, with the effect of $\mu_2$ being more pronounced. Figure 3(b) demonstrates that an increase in $\mu_2$ reduces $I_2$, while $\mu_1$ has a negligible impact. Figure 3(c) indicates that both $\mu_1$ and $\mu_2$ contribute to a reduction in the total number of infected individuals, although the effect of $\mu_2$ is more pronounced. The effect of the impaction of limited medical resources, $h_1$ and $h_2$, on the numbers of sensitive infected individuals, insensitive infected individuals, and total infected individuals at steady state is shown in Figure 4. As shown in Figure 4(a), a decrease in $h_1$ and an increase in $h_2$ both lead to a reduction in $I_1$. Figure 4(b) indicates that $I_2$ decreases as $h_2$ decreases, while $h_1$ has almost no impact. Figure 4(c) demonstrates that reductions in both $h_1$ and $h_2$ contribute to a decrease in the total number of infected individuals. In general, while the impact of different parameters varies, increasing the influence of media information and medical resources contributes to a reduction in the total number of infected individuals.

Next, we examine the impact of several key parameters on the total number of infections and model stability, as shown in Figure 5. It is evident from the figure that variations in these key parameters significantly influence both the changes in and the stability of the total number of infected individuals. Then, we analyze the relationship between changes in key parameters and stability through the following bifurcation diagrams.

Figures 6 and 7 are Hopf bifurcation diagrams of $I_1$ and $I_2$ caused by parameters. As shown in Figure 6(a), for the range of the reaction of sensitive susceptible individuals to the media, $0 < {\alpha _1} \le 0.0102$, the disease is asymptotically stable. For $0.0102 < {\alpha _1} \le 0.043$, the periodic oscillations (limit cycle) will appear, i.e., the disease outbreak will occur repeatedly. However, for a higher value of $\alpha_1$, i.e., for ${\alpha _1} > 0.043$, the periodic solutions disappear and the disease again becomes stable. For Figure 6(b), when $0 < {\alpha _2} \le 0.014$, the disease is asymptotically stable. When ${alpha _2} > 0.014$, periodic oscillation (limit cycle) will occur, that is, disease outbreaks will occur repeatedly. For Figure 6(c), when $0 < {\mu _1} < 1$, there will always be periodic oscillation (limit cycle), and the disease will break out repeatedly. In Figure 6(d), when $0 < {\mu _2} \le 0.16$, the disease will remain stable, while when $0.16 < {\mu _2} \le 1$, the disease will occur in periodic outbreaks. In Figure 6(e), when $0 < {h_1} < 4$, the disease will occur periodic outbreaks. The situation of the Figure 6(f) is similar to that of Figure 6(a). When $0 < {h_2} \le 0.006$, the disease is stable, while when $0.006 < {h_2} \le 0.068$, the disease is unstable and periodic oscillation occurs, and when $h_2$ is increased again, the disease will become stable.

Figure 7(a)–(f) analytical methods parallel to those employed for Figure 6. There are also several situations caused by parameter changes, such as from stable to unstable, unstable to stable, unstable in the whole stage, and from stable to unstable and then to stable. Therefore, media information and medical resources play crucial roles in inducing complex dynamics in model (2.2).

5.2.   Numerical simulation of bifurcations

Due to computational challenges, we use the MATCONT package to perform numerical simulations and analyze the dynamical behaviors and bifurcations of model (2.2) under various conditions. We shall illustrate forward bifurcation, backward bifurcation, and forward-backward bifurcation in the following.

A forward bifurcation is illustrated in Figure 8, and Figure 8(a) gives variation in equilibrium level of sensitive infected individuals with the basic reproduction number (associated with parameter $\Lambda _1$) and corresponding bifurcations, say, when $\mathcal{R}_0 = 1.0214$, the model (2.2) undergoes Hopf bifurcation. Figure 8(b) is a partial enlarged view of Figure 8(a) near $\mathcal{R}_0 = 1$. Figure 8(c)–(f) illustrates the time series of the model (2.2) with different parameters $\Lambda_1$ which stabilize to either the disease-free equilibrium or endemic state or bifurcated periodic solution. The first Lyapunov coefficients at both Hopf points are negative, which means that the Hopf bifurcation is supercritical and the periodic orbits are born stably, that is to say, the model exists a stable limit cycle when $1.0214 < \mathcal{R}_0 < 1.413$. Also, we can see that when $\mathcal{R}_0 < 1$, the model (2.2) only has the DFE, which is globally asymptotically stable, shown in Figure 8(c); when $1 < \mathcal{R}_0 < 1.0214$ and $\mathcal{R}_0 > 1.413$, the model (2.2) has an unstable DFE and a globally stable endemic equilibrium, shown in Figure 8(d),(f); when $1.0214 < \mathcal{R}_0 < 1.413$, the model (2.2) has an unstable DFE and an unstable endemic equilibrium, but has a bifurcated periodic solution around this unstable endemic equilibrium, shown in Figure 8(e).

A backward bifurcation is shown in Figure 9, and Figure 9(a) illustrates the saddle-node bifurcation at $\mathcal{R}_0 = 0.585$ and backward bifurcation at $\mathcal{R}_0 = 1$. Figure 9(b)–(d) are time serious diagrams corresponding to different $\mathcal{R}_0$ of backward bifurcation diagram. When $\mathcal{R}_0 < 0.585$, model (2.2) has a unique stable DFE, shown in Figure 9(b). As $\mathcal{R}_0$ increases, the $E_0$ coexists with a stable endemic state ${E^*}$ for $0.585 < \mathcal{R}_0 < 1$, shown Figure 9(c), and further increasing $\mathcal{R}_0$ leads to model (2.2) stabilizing to the stable endemic equilibrium for $\mathcal{R}_0 > 1$, shown in Figure 9(d).

The forward-backward bifurcation is illustrated in Figure 10(a) where saddle-node bifurcation occurs for $\mathcal{R}_0 > 1$ and Hopf bifurcation may also occur, and Figure 10(b) is a local magnification of Figure 10(a) around $\mathcal{R}_0 = 1$. We find two Hopf bifurcation points at $\mathcal{R}_0 = 0.6451$ and $\mathcal{R}_0 = 1.0092528$ and two saddle-node bifurcation points at $\mathcal{R}_0 = 0.617$ and $\mathcal{R}_0 = 1.0288$. The first Lyapunov coefficients at those two Hopf points are both positive, which means that the Hopf bifurcations are subcritical and the bifurcated periodic orbits are born unstable. To further analyze how periodic orbits change, we plot local bifurcation diagrams (Figure 10(c),(d)) and the variation in period of cycles versus $\mathcal{R}_0$ (Figure 10(e),(f)). Figure 10(c),(d) show that once the periodic orbit appears, the minimum value of $I(t)$ of the periodic orbit is infinitely close to the value of $I(t)$ at the unstable saddle point (or DFE) at $\mathcal{R}_0 = 0.6534$ and $\mathcal{R}_0 = 1.009208$. Figure 10(e),(f) shows that the periods eventually go to infinity at $\mathcal{R}_0 = 0.6534$ and $\mathcal{R}_0 = 1.009208$, at which the homoclinic bifurcations occur and the limit cycle is replaced by a homoclinic orbit.

In summary, the numerical bifurcation diagrams reveal the occurrence of forward bifurcation, Hopf bifurcation, backward bifurcation, forward-backward bifurcation, and homoclinic bifurcation. During the initial phase of an emerging infectious disease outbreak, intensifying public awareness campaigns and implementing containment measures can effectively reduce the initial number of infections. This reduction enables the system to stabilize the infection count at a lower endemic equilibrium, thereby mitigating long-term public health burdens.

5.3.   Constant control and optimal control with limited resources

To investigate the information allocation between two populations with different information reception capabilities, numerical simulations will be conducted to determine the total number of infections and the constant media resource allocation rate $u$ under the constraint of limited media resources. Subsequently, a time-varying optimal allocation rate $u(t)$ control strategy will be proposed to minimize both the total number of infections and the cost associated with media information allocation. In this section, we use the media resource data from ^[38] as the parameter values for $V$. In their study, the daily news data from Xinhuanet (news.cn) ranged approximately between 5 to 45 articles per day. Therefore, we also adopted this data range in our analysis, where 5 articles per day indicate insufficient media coverage (limited resources), and 45 articles per day represent sufficient media coverage (abundant resources).

First, Figure 11 illustrates the relationship between the constant allocation rate (u) of model (4.1) and the total number of infected individuals across varying resource levels, considering two groups with three different types of heterogeneity. From Figure 11(a), it can be observed that when heterogeneity is low and $V \le 25$, media resources should be allocated to the group with higher information sensitivity. However, when $V > 25$, resources should be allocated to the insensitive group. Figure 11(b) indicates that when the heterogeneity between the two groups is moderate, the allocation strategy shifts at $V = 25$, with resources now being allocated to the insensitive group. In cases of high heterogeneity between the two groups, as shown in Figure 11(c), media resources should also be allocated to the insensitive group when $V = 15$. Therefore, the allocation of media resources is closely related to the heterogeneity between the two groups. As heterogeneity increases, resources are more frequently allocated to the insensitive group.

Next, to gain a deeper understanding of the optimal media allocation strategy for minimizing infections and associated costs, we simulated the optimal control in model (4.1) and derived the corresponding optimal control. In the context of three different levels of group heterogeneity, we use the same parameter values, except for the values of $\alpha_1$ and $\alpha_2$, and apply the forward-backward sweep method to solve the control problem. We numerically calculated optimal control strategies based on the iterative method used in ^[39]. This algorithm addresses the optimal control problem by simulating the system dynamics through Ordinary Differential Equation Solvers and updating the control strategy using adjoint equations. During each iteration, the objective is to minimize the cost function by adjusting the control variables. After every iteration, the control variables are updated and the system states are recalculated until the convergence criteria are satisfied. The specific steps are as follows:

1) Parameter initialization:

● Set global variables such as $\lambda_1$, $\lambda_2$, $\beta_{11}$, $\alpha_1$, etc., which are parameters related to the model.

● Initialize control variables, iteration parameters, such as $u$, $u{(t)}$ (current and previous iteration control variables), and initial system states $y_0$.

2) Model definition:

● Define different models corresponding to different control strategies.

● These models describe the system state evolution over time with control inputs affecting the system dynamics.

3) Objective function definition:

● The objective function $J$ is used to evaluate the effectiveness of the current control strategy. It is computed by summing the cost at each time point, which usually involves factors such as infection rates and control cost.

● During each iteration, calculate the new objective function and compare it with the previous iterations value to check for convergence.

4) Control input calculation:

● Update control input $u$ based on the adjoint equation (Lagrange multiplier) and the current state $y$.

● Control input $u$ is updated using the adjoint variable ${\lambda _{{S_1}}}, {\lambda _{{I_1}}}, {\lambda _{{R_1}}}, {\lambda _{{S_2}}}, {\lambda _{{I_2}}}, {\lambda _{{R_2}}}, {\lambda _M}$ and the system state $y$.

● Ensure the control input is within a specified range, and update the control variable using a step size.

5) Adjoint equation solution:

● Solve the adjoint equation to update the adjoint variables ${\lambda _{{S_1}}}, {\lambda _{{I_1}}}, {\lambda _{{R_1}}}, {\lambda _{{S_2}}}, {\lambda _{{I_2}}}, {\lambda _{{R_2}}}, {\lambda _M}$. This step is necessary in optimal control problems.

● In the backward step, starting from the terminal time, the adjoint equations (typically linked to the Lagrange multiplier equations of the optimization problem) are solved to update the adjoint variables.

6) Optimization iteration:

● Iterate multiple times, updating the control variable $u$ at each step until the objective function changes by less than a threshold or the maximum number of iterations is reached.

● After each iteration, print the current control variables $u$ and state variables $y$, and update the objective function values.

7) Visualization:

● During the iterations, plot various graphs, such as the infection rate change and the evolution of control strategies.

● Display the results for different control strategies.

As shown in Figure 12(a),(c),(e), the optimal control strategy is implemented from the onset of the epidemic. Initially, media resources are gradually allocated to sensitive groups. As the disease progresses, all media resources must be allocated to sensitive groups. However, as the infection rate within the population declines, media resources should gradually be shifted toward less sensitive groups. The key distinction across varying levels of population heterogeneity is that, as heterogeneity increases, the timing of the transition in media resource allocation occurs progressively earlier.

Figure 12(b),(d),(f) presents the infected humans prevalence, where the red dashed line represents the infection rate with a control strategy of 0.2, the blue dashed line represents the infection rate with a control strategy of 0.5, the purple dashed line represents the infection rate with a control strategy of 0.8, and the green solid line indicates the infection rate under the optimal control strategy. It is evident that under optimal control, the infected humans prevalence is minimized. However, the overall difference in the effectiveness of infection prevalence control across the three types of heterogeneous populations is minimal. Only in groups with high heterogeneity does controlling the infected humans prevalence below $20\%$ lead to slightly more effective control compared to the other two strategies. Thus, the distribution of media resources follows a heterogeneity-dependent allocation paradigm: under low heterogeneity conditions, priority is given to information-responsive groups, whereas higher heterogeneity triggers a proportional reallocation to demographic segments exhibiting limited responsiveness to media interventions.

6.   Conclusions and discussion

In this study, we proposed an improved infectious disease $SIR$ model that divides the population into two distinct groups: the information-sensitive group and the information-insensitive group. This model incorporates two nonlinear functions to examine the impact of media coverage and limited medical resources on disease transmission. It more accurately reflects the actual dynamics of infectious disease transmission under varying heterogeneity conditions. We analyze the well-posedness of the original model, but due to its high nonlinearity, we focused on the case where only media information influences and investigated the conditions for disease extinction and persistence. Numerical simulations were conducted on media and medical resource parameters. The bifurcation analysis revealed various dynamic behaviors, including forward bifurcation, backward bifurcation, Hopf bifurcation, saddle-node bifurcation, and homoclinic bifurcation. In the control problem, we further considered the allocation of limited media resources and explored the issue from both constant and optimal control perspectives. Under constant control, we simulated the changes in the number of infected individuals for different types of population heterogeneity, considering varying media allocations and limited media resources. In the case of optimal control, we applied optimal control theory to explore the problem related to information allocation and derived optimal control strategy. In addition, we validated the effectiveness of the proposed control strategy by simulating the optimal media resource allocation for different types of population heterogeneity, and compared the similarities and differences in control strategies and effects among different types of groups.

Our model incorporates factors such as media coverage and medical resource limitations, particularly in the early stages of a pandemic, where these factors significantly influence disease transmission. It also captures the variations in media responses across different populations, which is essential for developing targeted public health policies. By analyzing simplified conditions for disease extinction and persistence, we provide theoretical insights into epidemic dynamics. However, media coverage and medical resource limitations do not affect the threshold conditions for disease extinction or persistence. Numerical simulations of media and medical resource parameters illustrate their impact on the number of infected individuals. Overall, increasing public responsiveness to media and medical resources helps reduce infections, and both factors influence the stability. Bifurcation analysis reveals complex dynamics under varying parameters, highlighting the challenges of controlling disease spread. Two control methods were simulated with limited media resources. Under constant control, resources tend to be allocated to insensitive groups as population heterogeneity increases. In optimal control, a key distinction across varying levels of heterogeneity is that as heterogeneity increases, resources gradually shift toward insensitive groups.

Compared to ^[21] and ^[26], our work extends the model by jointly considering the effects of media coverage and medical resources, and explores control strategies for populations with varying heterogeneity. Our research suggests that during the early stages of an emerging infectious disease outbreak, public health organizations can reduce the initial number of infections and stabilize the infected population at a lower level through measures such as quarantine controls or intensified disease awareness campaigns. Furthermore, the allocation of media resources is closely tied to population heterogeneity: when heterogeneity is low, resources are predominantly allocated to information-sensitive groups, while as heterogeneity increases, resources gradually shift toward insensitive groups. Thus, our study may contribute to advancing their work in this area. However, certain limitations remain, particularly in the theoretical analysis, where the high nonlinearity of the model constrained us to focus on the extinction and persistence conditions of a simplified model. Additionally, our current framework does not account for the heterogeneity of media information, such as its veracity (truthfulness vs. misinformation) and practical relevance (actionable guidance vs. speculative content), which are critical factors requiring explicit consideration in real-world public health communication strategies. This aspect can be further refined in future research.

Use of AI tools declaration

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

Acknowledgments

The authors are very grateful to the editor and the anonymous reviewers for their careful review and valuable advice, which helped to improve the presentation of the paper. The authors sincerely thank Dr Xiulei Jin, Ning Wang, and Wei You for their valuable contributions to the numerical simulation section of this paper, which significantly improved the quality of this work. The authors acknowledge support from the National Natural Science Foundation of China (No.12271401, W2421101) and the Natural Science Foundation of Tianjin, China (No.22JCYBJC00080).

Conflict of interest

The authors declare there is no conflict of interest. Shengqiang Liu is a guest editor for Mathematical Biosciences & Engineering and was not involved in the editorial review or the decision to publish this article. All authors declare that there are no competing interests.