Effects of Media Coverage on Global Stability Analysis and Optimal Control of an Age-Structured Epidemic Model with Multi-Staged Progression

: In this paper, a hybrid SEIAM model for infectious disease with a continuous age structure is established, where combined dynamic effects of media coverage and multi-staged infected progression on threshold dynamics are discussed. Sufﬁcient conditions for uniform persistence of the solution are studied by using the basic reproduction number. By constructing appropriate Lyapunov functions, the global stability analysis of endemic equilibrium is investigated based on Lyapunov–LaSalle’s stability theorem. In order to minimize costs incurred due to applied controls and infectious disease burden, an optimal cost-effective control strategy associated with disease treatment and media coverage is discussed. Numerical simulations are carried out to show consistency with theoretical analysis.


Introduction
According to some recent epidemiological data collected in medical observations, it is found that parasite levels in vectors or viral loads in an infected population have a great influence on disease infectivity [1], and not all of the infected population behave in a homogeneous way during the entire infectivity period [2][3][4]. Some infectious diseases with a comparatively long infectivity period, such as HIV [1], bilharzia [2], HBV and HCV [3,4], usually exhibit various infecting forces at successive discrete stages of infection. Some staged progression models are proposed to depict disease progression with typical varying characters in infectivity transmission [5][6][7]. For epidemic models with multiple infectious stages and general distributions for the stage durations (exponential distribution [5] and Gamma distribution [6] as special cases), it is shown that disease dynamics are determined by the derived basic reproduction number. The global threshold dynamics and uniform persistence of stage-processed disease models are discussed based on the basic reproduction number in [7].
In recent decades, many research efforts have been made to discuss the potential impact of age on disease progression, and some age-structured models have also been established to depict the vaccination or latency phenomena in the studied epidemiology [8,9]. Although these mathematical models are usually composed of partial differential equations and corresponding dynamical analysis is particularly challenging, global stability for some type of age-structured models have been investigated by constructing Lyapunov functionals in [10][11][12] and the references therein. Complex bifurcation phenomenon stability analyses are discussed in some multi-group epidemic models with age structure [13][14][15][16] and the references therein. A novel multi-group SVEIAR epidemic model with ages of vaccination and latency is formulated in [17], where the sufficient conditions for the occurrence of backward bifurcation in epidemic models are investigated. In order to discuss the possible effects of the variable infected population in each stage of disease dynamics, some epidemic model multi-age structures are proposed in [18][19][20], and the uniform persistence of the models and global stability analysis around the epidemic equilibrium are investigated. In [21], the authors propose a reaction-diffusion SIR epidemic model with multi-group and the nonlinear infectious force function is investigated, and global asymptotic stability for disease-free and endemic equilibria is discussed by using the Lyapunov functions method.
It is well known that media coverage reduces the prevalence of infectious disease and shortens the duration of infection [22]; public health alerts through social media and TV advertisements play an important role in effectively communicating with people regarding the spread of any infectious disease during a short period of time [23]. The impact of public awareness and local behavioral response is studied in [24], and the awareness among the infected population reduces the infectivity due to pharmaceutical interventions. Since human behavior is usually changed toward the diseases through mass media, it is constructive to prevent the spread of infectious diseases as the public relies on what mass media projects to them [25,26]. By assuming the public health alert programs are implemented proportional to the infected population, there has been recent progress on cost-effective analysis for the social media coverage to promote vaccination against infectious disease [27,28], and some recent studies were conducted to discuss the increment in vaccination coverage due to social media and TV advertisements, such as [29][30][31][32], by considering media alerts as a decreasing function of the infected population.
Recent theoretical research shows that vaccination and treatment are the most effective as well as applicable preventive and control strategies during the outbreak of infectious disease. In order to evaluate optimal strategies of social distancing and vaccination, an age-structured dynamic system is established to describe the transmission mechanism of seasonal influenza in [33]. By using Pontryagins' maximum principle, the authors identify the best way of reducing morbidity and mortality at a minimal cost. By choosing the treatment (a pharmaceutical intervention) and constructive public health alerts (a nonpharmaceutical intervention) as control tools, an optimal control problem is formulated in [34], which is utilized to minimize the cost and disease fatality. In [35,36], the author discusses the optimal distribution strategies of an infectious disease with the mixing of two sub-groups under a resource-constrained environment. An optimal control strategy is formulated [37] to explore cost-effective intervention strategies associated with vaccination and water sanitation, which are designed for the prevention and intervention of cholera epidemics. By considering information-induced vaccination and treatment as controls, optimal control profiles for the applied controls are obtained in [38]. For some specific infectious diseases (foot-and-mouth disease model in [39] and liver cirrhotic disease model [40]), optimal vaccination and treatment measures are formulated.
To the author's best knowledge, there are few results concerning mathematical modeling of the age-structured epidemic model with stage progression and media coverage. The combined dynamic effects of multi-staged infected progression and public health alerts on epidemiological dynamics have not been investigated before. The optimal control strategy associated with treatment and media coverage has not been discussed in the related references mentioned above.

Model Formulation
Based on the above analysis, some hypotheses are proposed as follows: Hypothesis 1 (H1). The total population at any time t is divided into the susceptible population S(t), exposed population E(t), time and age-dependent infected population I(t, a) and aware population A(t). It is assumed that media campaigns induce behavioral changes among the susceptible population, which form an isolated aware class fully protected from the infectious disease. Λ denotes the constant recruitment rate of the susceptible population. d 1 , d 2 and d 4 denote the per capita death rate of S(t), E(t) and A(t), respectively. δ stands for per capita transfer rate from the exposed population to the infected population, c denotes the dissemination rate of awareness among the susceptible population, h 0 represents the transmission rate of the aware to the susceptible population.
Hypothesis 2 (H2). In order to investigate the possible effects of multiple infectious stages and variable infectivities in each stage on the epidemiological dynamics, the infection age in the k-th stage is assumed to be characterized by its relationship to the corresponding adjacent stages. The entire infectious period [0, +∞) in this paper is assumed to be partitioned into (k + 1)-th stages defined by infection age intervals (a k−1 , a k ] (k = 1, 2, · · · , n) such that 0 = a 0 < a 1 < · · · < a n < a n+1 = ∞.

Hypothesis 3 (H3).
It is assumed that the susceptible population is vulnerable to the infectious disease as well as public health alerts regarding social media advertisements and TV. µ k (a) represents the age-dependent transmission coefficient in the k-th infectious stage, which describes the varying probability of infectiousness in the k-th stage. ν k (a) denotes the age-dependent recovery rate. d 3k stands for the per capita death rate of the infectious population at the k-th stage, β k (a) denotes the transmission rate of the infectious population from the k-th stage to the next (k + 1)-th stage.
Hypothesis 4 (H4). M(t) represents the cumulative number of social media and TV advertisements at time t. It is plausible to assume that the cumulative number of social media and TV advertisements proportionally increases with the number of infected individuals. It is also assumed that the growth rate of social media and TV advertisements is proportional to the number of the infected population with a decreasing function of the aware population. l represents the half-saturation point for the impact of social media and TV advertisements on the susceptible population. M 0 and η represent the baseline number and diminution rate of social media and TV advertisements, respectively. r k (a) denotes the age-dependent growth rate of social media and TV advertisements. θ k (a) denotes the age-dependent decay rate in social media and TV advertisements due to an increase in the aware population. ω k (a) represents the half-saturation point as r k (a)θ k (a) A(t) ω k (a)+A(t) attains half of its maximum value r k (a)θ k (a). Hypothesis 5 (H5). All constant parameters are nonnegative, and the age-dependent functions are all nonnegative, bounded and integrable in their definition intervals.
Keeping all these Hypotheses (H1)-(H5) in mind, we formulate a multi-stage SEIAM model for infectious diseases with media coverage and multi-staged infection progression, which is governed by the following hybrid model with ODE and PDE, for k = 1, · · · , n + 1 with a 0 = 0, a n+1 = ∞ and β n+1 (a) = 0 for any a ∈ (a n , a n+1 ). The boundary conditions are For any t > 0, k = 1, · · · , n + 1 and initial conditions

Remark 1.
Recently, some epidemic models with media-induced nonlinear incidence and agedependent vaccination have been established in [41,42]. Uniform persistence and the sharp threshold dynamics are investigated, and the efficacy of vaccines is discussed due to the variation in vaccination reproduction number. However, the combined dynamic effects of multiple disease progressions and public health alerts on epidemiological dynamics have not been investigated. Compared with the work in [41,42], model (1) can not only describe infectious disease progression through multiple stages but also assess the dynamic impact of social media and TV advertisements on the spread of an infectious disease. The treatment measures generally minimize the serious consequences of infection, and efficient media coverage act as a constructive public health alert, which reduces the infection rate. With these proper control measures, the authorities may significantly reduce wealth loss and mortality due to disease burden. Compared with recent references [33][34][35][36][37][38][39][40] concerning the optimal control strategies for infectious disease, the treatment and diminution rate of media coverage are first adopted as control variables in this paper. An optimal cost-effective control strategy is discussed to minimize the costs incurred due to applied controls and infectious disease burden.
The rest of the sections of this paper are organized as follows. In the Section 2, nonnegativity and the uniform persistence of the solution are studied by using the basic reproduction number. By constructing the appropriate Lyapunov functions, the global stability analysis of endemic equilibrium is investigated. In the Section 3, an optimal costeffective control strategy is discussed to minimize costs incurred due to applied controls and infectious disease burden. In the Section 4, numerical simulations are provided to support the theoretical analysis. Finally, this paper ends with a conclusion.

Qualitative Analysis of SEIAM Model
In this section, the combined dynamic effects of media coverage and multi-staged infected progression on infectious disease dynamics are discussed. Sufficient conditions for uniform persistence of the solution are studied by using the basic reproduction number. By constructing appropriate Lyapunov functions, global stability analysis of the endemic equilibrium is investigated based on the Lyapunov-LaSalle stability theorem.

Basic Reproduction Number and Endemic Equilibrium
In order to describe the age-specific survival probability of an infected person in the age interval (a k−1 , a k ], we define and some mathematical notations are also defined for k = 1, · · · , n + 1, Remark 2. It should be noted that V k (k = 1, 2, · · · , n + 1) represents the probability that an infected individual survives during the k-th infection age and proceeds to the (k + 1)-th infection age. (δ+d 2 )[cM 0 +(h 0 +d 4 )(l+M 0 )] is the fraction of individuals surviving through the exposed stage and reaching the first infection age. Q k Π k−1 j=1 V j stands for the part of secondary infections generated from the k−th infection age, when a typical infected individual contacts a susceptible individual. Hence, the sum of n + 1 terms ∑ n+1 k=1 Q k Π k−1 j=1 V j denotes secondary infections from all n + 1 infection ages.
By using Remark 2 and simple computations, it is easy to show the basic reproduction number is as follows: Based on system (1), if R 0 > 1, then it follows from simple computations that the endemic equilibrium P * = (S * , E * , I * 1 (a), · · · , I * n+1 (a), A * , M * ) is as follows,

Nonnegativity and Uniform Persistence of Solution
In order to facilitate the following proof, vector w is defined as follows, it is a solution vector, which is endowed with the following norm, Proof. Based on the standard theory for age-structured models [43], boundary conditions (2) and initial conditions (3), by integrating I k (a, t) on (a k−1 , a k ] it yields, According to the classical existence and uniqueness results for functional differential equations, by substituting (8) into system (1), it is easy to show that I k (a, t) is nonnegative for any nonnegative initial value, and there exists a unique solution for system (1).
If there exists t 1 such that S(t 1 ) = 0 and S(t) > 0 holds for 0 < t < t 1 , then it implies that dS(t 1 ) dt > 0 and S(t) ≥ 0 holds for t ≥ 0. By using similar arguments, the nonnegativity of other variables, E(t), R(t), A(t) and M(t), can be shown for any nonnegative initial values. This completes the proof.
In order to guarantee the existence of the semiflow solution for system (1), which admits a compact global attractor, system (1) can be rewritten as the following Cauchy problem defined on boundary condition (2), where A : Dom(A) → R 4 × L 1 (0, +∞), R n+1 × R n+1 , z = (w, 0, · · · , 0) T ∈ R 4 × L 1 (0, +∞), R n+1 × {0} n+1 take the following form, with z(0) = z 0 ∈ B + , Dom(A) = z ∈ B + u k (·) ∈ Z 1,1 [0, +∞), k = 1, · · · , n + 1 , Z 1,1 represents a Sobolev space and B + is defined as follows, Furthermore, F is defined as follows, . . a n a n−1 β n (a)u n (a, t)da It follows from Theorem 2.1 and Theorem 2.2 in [44] that there exists a unique solution semiflow Z(t) : B + → B + defined by system (9). Let , it follows from simple computations and system (1) that By using (13) and the comparison principle [45], it can be obtained that which reveals that all solutions of system (1) and the corresponding system (9) are uniformly bounded. Hence, the solution semiflow Z(t) is point dissipative. According to Theorem 3.6.1 [45], it can be obtained that Z(t) is compact for all t > 0, and it is easy to show that Z(t) has a compact global attractor in B + based on Theorem 1.1.3 [46]. Based on the above analysis, Lemma 2 can be concluded as follows. When t is large enough, Z(t)v →v holds for each x ∈ ∂I, where I, ∂I andv are defined as follows, Proof. Based on (14), it is easy to show that S(t) ≤ Λ+η M 0 α , for any initial point in ∂I and   S 0 , I 0 1 , · · · , I 0 n+1 , 0, · · · , 0 n+1 , 0, 0   ∂I, we consider the following auxiliary system By using the comparison principle [46], it can be obtained that is a solution of system (15).
Furthermore, there exists a compact subsect of I 0 ⊂ I, which is also a global attractor for {Z(t)} t≥0 in I.
Proof. According to Lemma 3, it is easy to show thatv is globally asymptotically stable in ∂I. In order to investigate the asymptotic behavior of the solution to system (15) from the neighborhood ofv, wherev is defined in Lemma 3.
For the following auxiliary system (20) with k = 1, · · · , n + 1 and J ≥ 0, According to the comparison principle [46], it yields that where E J (t), (22), By similar arguments in [18,47], the following notations are defined to facilitate discussing the characteristic equation of system (22), where a ∈ (a k−1 , a k ], k = 1, · · · , n + 1. When R 0 > 1 and there is a sufficiently large J > 0, it follows from simple computations that the dominant eigenvalues of system (22) satisfy the following characteristic equation It follows from (23) (22) can be linearized as follows: where , and the closed linear operator B is as follows, Rev . According to the Hille-Yosida theorem [48], a strongly continuous semigroup { Z(t)} t≥0 can be generated by the operator (B, D(B)). For any v J ∈ R 4 + × I × {0} n+1 from the neighborhood ofv, Φ J denotes the specific projector on the eigenspace with respect to the dominant eigenvalue v * J , it follows from the similar arguments in [18,49], it yields According to the recurrence relationships of ϕ(t), it yields that and By substituting (28) into (27) and (29), ϕ 1 , ϕ n+3 and ϕ n+4 can be obtained as follows, dv J > 0 and I 0 1,J , · · · , I 0 n+1,J ∈ I, it follows from (28) and (30) Hence, it follows from the above analysis that (31) contradicts (28). Based on Theorem 4.2 in [50], it is easy to show that the semiflow {Z(t)} t≥0 is uniformly persistent with respect to (R 4 + × I × {0} n+1 , ∂I ), i.e., there is a ζ > 0 such that there is a solution of system (9) with an initial value (S 0 , E 0 , I 0 1 , · · · , I 0 n+1 , 0, · · · , 0, A 0 , M 0 ) ∈ R 4 + × I × {0} n+1 , and Furthermore, according to the above analysis, it is easy to show that there exists a compact subsect of I 0 ⊂ I, which is also a global attractor for {Z(t)} t≥0 in I.

Stability Analysis of Endemic Equilibrium
Theorem 2. When R 0 > 1, the endemic equilibrium P * is locally stable.
Proof. According to (7), there exists an endemic equilibrium P * when R 0 > 1, by defining some following new variables with exponential form, and By substituting (32) into system (1), it yields that ∂I k (a, t) ∂a By separating and solving Equation (34), I * k (a) andĪ(a) are obtained as follows, where g k (a) and U k are defined in (5), andg k (v, a) is defined in (23). By substituting (33) into system (1),S,Ē,Ā andM are obtained as follows, According toĪ 1 (a 0 ) = δĒ, the characteristic equation is expressed as follows: In order to investigate the possibility of a solution to Equation (37) having a nonnegative real part, i.e., Re(v) ≥ 0, some discussions are provided as follows.
If v = 0, then Equation (37) is rewritten as which follows that Furthermore, if Re(v) > 0, then it is easy to show which also leads to a contradiction. Hence, it can be concluded that Equation (37) only has eigenvalues with negative real parts and the endemic equilibrium P * is locally stable. Theorem 3. When R 0 > 1, if C 1 > C 2 and C 3 > 0, then the endemic equilibrium P * is globally asymptotically stable.

Optimal Control Strategy
In this section, we formulate and analyze an optimal control strategy considering the diminution of media coverage and treatment to control the disease progression and balance the economic loss generated due to infectious disease and treatment measures. We will adopt time and age-dependent treatment ρ(a, t) on the infected population and a time-dependent diminution rate of media coverage η(t) on the cumulative number of social media and TV advertisements. Hence, the controlled system takes the following form, where system (49) shares the same boundary and initial values with system (1), other parameters share the same mathematical forms and biological interpretations described in system (1). The proper treatment measures and efficient media coverage significantly minimize the serious consequences of infection associated with wealth loss and mortality due to disease burden. Our aim is to optimize the use and minimize the cost of these two mentioned measures, and the optimal control problem is formulated as follows: where G and H(ρ(a, t), η(t)) are defined as follows, where ρ max and η max are positive upper constants of ρ(a, t) and η(t), respectively, and where H(ρ(a, t), η(t)) = T 0 n+1 ∑ k=1 a k a k−1 γ 11 I k (a, t) + γ 12 ρ(a, t)I k (a, t) + 1 2 γ 21 ρ 2 (a, t)da dt (52) where state variables are all subject to system (49), γ 11 and γ 13 are positive weights with respect to disease infection and the cumulative amount of media coverage. γ 12 represents the positive unit cost of controls at different levels. γ 21 and γ 22 denote positive constants associated with two quadratic terms accounting for control costs of disease infection and media coverage. In order to derive the optimal problem (50) with respect to time-dependent and age-dependent control functions, some derivations of the adjoint system (53) and the characterization of the optimal control functions are discussed as follows.
where σ denotes a scalar, φ is a constant balancing of the different units between time and age and f 1 (a, t) and f 2 (t) are arbitrary functions. ρ(a, t) and η(t) in system (49) are replaced with ρ(a, t) + σ f 1 (a, t) and η(t) + σ f 2 (t) in adjoint system (53), respectively. Since the solutions of system (53) are associated with σ, we will differentiate system (53) with respect to σ at σ = 0 and the system for sensitivities are obtained as follows, where F S , F E , F I k , F A and F M denote the derivatives of S, E, I k , A and M with respect to σ and are evaluated at σ = 0 (k = 1, · · · , n + 1), respectively.
It is easy to show that adjoint system (60) has a mild solution with Lipschitz continuity, Theorem 4. If (ρ * , η * ) is a pair of optimal control subjects to optimal control problem (50), then the characterization of optimal control (ρ * , η * ) takes the following form Based on (62), the characterization of the optimal control is as follows, Consequently, it follows from the above analysis that a pair of optimal controls can be obtained as follows: ρ * = min{max{0,ρ}, ρ max }, η * = min{max{0,η}, η max }.

Remark 3.
Generally, it is difficult to solve the sequence of optimal controls and the associated state variables converging with the optimal controls due to a lack of compact characters in the age-structured first-order PDEs. Hence, Ekeland's principle [51] will be utilized to minimize the sequences of the approximate functionals.
According to the corresponding analysis in Chapter 14 [52] that there exists a pair of controls , and the optimal control objective functional is defined as follows, If (ρ ε , η ε ) is a pair of optimal control subjects to the optimal control problem minimizing H ε (ρ, η), then it follows from the corresponding analysis in Chapter 14 [52], it is easy to show k=1 (a k−1 , a k ] × (0, T).
is sufficiently small, then there exists a pair of unique optimal controls (ρ * , η * ) subject to optimal control problem (50).
Proof. First, two functions are defined as follows, , ρ max , For two pairs of controls, (ρ 1 , η 1 ) and (ρ 2 , η 2 ), according to L ∞ Lipschitz estimations for the states and adjoint functions, it yields that where G i (i = 1, 2) are positive Lipschitz constants associated with the essential boundedness of the state variables and adjoint functions.

Numerical Simulation
In this section, some numerical simulations are conducted to show the feasibility of our analysis in the previous section regarding global stability and optimal control strategy, where values of parameters in numerical simulations are partially taken from [18,31,37].
When progressing to the deterioration of diseases, it is well known that infectivity during the acute state is greater than that of the chronic one, and infection awareness among susceptible individuals increases due to obvious symptoms of chronic patients. In order to indicate varying degrees of transmissibility at the acute stage (from 0 to 2 months) and chronic stage (from 2 to 10 months), the age-dependent coefficients at different stages are as follows: 3.8e −0.7877a , 0 < a ≤ 2, 0, 2 < a ≤ 10. ν 2 (a) = 0, 0 < a ≤ 2, 0.025, 2 < a ≤ 10.
By using the set of parameter values, Λ = 1.7, d 1 = 0.025, d 2 = 0.01, δ = 1, l = 1200, h 0 = 0.025, η = 0.025, the basic reproduction number R 0 due to variations of c and M 0 are plotted in Figure 1a,b, respectively. From this figure, it may be noted that R 0 decreases with the increasing dissemination rate of awareness (c), as well as the increasing baseline number of social media and TV advertisements (M 0 ). R 0 becomes less than unity with the appropriate c and M 0 , which shows the importance of public health alerts in controlling the spread of infectious disease.
Dynamic impacts of infection age on the two successive infectious states (acute stage from 0 to 2 months, chronic stage from 2 to 10 months) are shown in Figures 3 and 4. The distributions of the infected individuals with respect to both age and time are plotted in Figure 3, and the corresponding dynamical responses of I 1 and I 2 with respect to age are plotted in Figure 4, which reflects the transmission heterogeneity of the infected individuals in acute (Figure 4a) and chronic infectious stages (Figure 4b).
In the following part, numerical simulations are carried out to show the effects of optimal controls on reducing infectious disease and minimizing the costs of control strategies, which are formulated and analyzed in the third section of this paper. By extending the forward-backward sweep method for ODE models, an iterative process is designed for numerical simulations of state and adjoint systems in the age-structured PDE model (60). In order to maintain relatively good numerical algorithm stability, a simple finite difference scheme is selected to facilitate the numerical implementations. By fixing the weight parameter values with respect to disease infection and cumulative number of media coverage as follows γ 11 = 3.5, γ 13 = 5, and varying the control costs by discussing three scenarios as follows, (i) low optimal cost controls: γ 12 = γ 21 = γ 22 = 0.35; (ii) moderate optimal cost controls: γ 12 = γ 21 = γ 22 = 3.5; (iii) high optimal cost controls: γ 12 = γ 21 = γ 22 = 35. The total number of infected individuals with respect to time at the acute stage is plotted in Figure 5a, and at the chronic stage is plotted in Figure 5b. The total number of infected individuals with respect to age at the acute stage is plotted in Figure 6a, and at the chronic stage is plotted in Figure 6b. It follows from the curves of Figures 5 and 6 that the infected individual with high control costs decreases quickly at the acute stage, which implies the spread of infectious disease can be significantly reduced with the increase in control costs. The infected individual with high control costs finally approaches a relatively steady state (approximately 0) at the chronic stage.

Conclusions
It is observed that social media and TV advertisements regarding the spread of infectious diseases generally have potential effects on behavioral changes among susceptible people as well as to control the spread of diseases. The infected population usually shows different infecting forces during various stages of progress. Hence, it seems more reasonable for the specific purpose of disease modeling to introduce infection age structure. In this paper, we propose a mathematical model to discuss how media coverage impacts the global stability and threshold dynamics of an infectious disease. In order to depict varying degrees of transmissibility at different stages of infection, the continuous age structure for each successive infectious stage during a long infective period is taken into account. Sufficient conditions associated with the basic reproduction number are derived for uniform persistence of solution in Theorem 1. By constructing appropriate Lyapunov functions, it follows from the Lyapunov-LaSalle stability theorem that global stability analysis of endemic equilibrium is investigated in Theorem 3.
As public health alerts spread among the susceptible population, behavioral changes such as protective measures and vaccination will help healthy individuals to avoid contracting infection. Consequently, we quantify the information-induced vaccination function as a control measure to vaccinate susceptible individuals so that the maximum number of individuals get vaccinated with minimum cost involved. Since the treatment includes diagnosis, medication, hospitalization and other related hygienic efforts depending on limited medical resources, the treatment function is chosen as another control for optimal treatment. The total incurred cost is determined by taking a weighted sum of the cost of productivity loss due to disease and costs incurred in applying control interventions. Analytical characterization and the existence of optimal control paths is carried out in Theorems 4 and 5, respectively. Author Contributions: Conceptualization, C.L.; methodology, C.L. and P.C.; software, Q.J.; validation, C.L., P.C. and L.C.; investigation, C.L., P.C. and Q.J.; resources, Q.J.; writing-original draft preparation, C.L. and P.C.; writing-review and editing, Q.J. and L.C.; supervision, L.C.; funding acquisition, C.L. and P.C. All authors have read and agreed to the published version of the manuscript.