Optimal control and stability analysis of an age-structured SEIRV model with imperfect vaccination

1.   Introduction

Imperfect or leaky vaccines are those that do not provide complete immunity to the disease in their hosts and are incapable of preventing transmission to others. It is known that leaky vaccines can reduce the symptoms of the disease but do not offer full protection against infection or the spread of the pathogen. Therefore, studying the impact of imperfect vaccines on the human population becomes crucial. Additionally, for certain infectious diseases, it has been observed that the effectiveness of the vaccine may diminish over time. In the recent COVID-19 pandemic, it has been noted that most vaccines remain effective for a specific duration before individuals become susceptible to the virus once again.

During the course of the follow-up period of six months, Stephen et al. ^[1] studied the vaccine efficacy of different vaccines (shown in Figure 1). Menni et al. ^[2] showed that in the case of COVID-19, the vaccine's efficacy remained high after five months among those under the age of 55. They also showed that the booster dose can restore the effectiveness of the vaccine. Jake et al. ^[3] also studied the vaccine efficacy and demonstrated that the evidence of waning immunity is overstated. They demonstrated that antibody levels can be boosted in the general population and robust clinical data are required to assess the need for additional doses. This naturally motivates us to consider a compartmental model for infectious diseases with a vaccine-waning rate.

Kermack and Mckendrick ^[5,6,7] developed ordinary differential equations (ODEs) based models to study the transmission of infectious diseases. These models are formulated based on the assumption that all individuals have identical vital rates regardless of age or size. Enough literature is available on ODE-based infectious disease models; more details can be found in ^{[8,9,10,11,12]}. In epidemic models, the "disease clock" is important for studying the long-term dynamics of an infectious disease. Moreover, in transmitting infectious diseases, the age of infection or the time since infection starts plays a vital role. Age-structured population models characterize individuals based on their age, which can be class-age or demographic age. Class age means the time since infection, and demographic age means the chronological age of individuals.

Iannelli ^[13] discussed age-structured population models and their role in describing the growth and interaction of populations. Thieme and Castilo ^[14] examined the role of infection age in the HIV viral dynamics. More details about age and size-structured population models and their applications in biology and epidemiology are given in ^[16]. For the first time, Arino et al. ^[17] discussed the role of leaky vaccines with a general waning function. They consider three classes S, I, and V (without age structure), corresponding to susceptible, infected, and vaccinated individuals. They showed that subthreshold equilibria are possible, which may play an important role while designing vaccine strategies. Optimal control strategies for age-structured compartmental models are discussed by ^[18,19]. Semi-random epidemic network models can also describe the dynamics of many infectious diseases ^[20]. Awel et al. ^[19] incorporated imperfect vaccination in their model for Malaria disease, and Wei et al. ^[21] studied imperfect vaccination on scale-free networks. Many authors have also considered delay terms in fractional-order population models with or without structure variables. For more details, we refer to ^{[22,23,24,25]}. There are many interesting works on age and size-structured models for infectious diseases (see, for example, ^{[26,27,28,29,30,31,32]}).

We propose an age-structured SEIRV model with imperfect vaccination and vaccine-waning rates. After formulation of our model, we show the existence and uniqueness of the solution. We use semigroup of operators to show the existence and uniqueness of the solution. We find steady-state solutions to our model and check the stability of the disease-free steady state. We also study the optimal control problem with vaccination rate as a control variable. Our objective is to minimize the cost of vaccination and to reduce the number of infected individuals. Using adjoint equations, we derive the form of the optimal control variable. We show that the optimal control variable depends on the vaccine waning rate, infection rate, and the cost of vaccination. We also solve our model numerically to validate our theoretical findings. Our study is mathematical, and we study the qualitative properties of the solution. We obtain an explicit form for the optimal control variable, which involves important parameters of the proposed model. Due to the complex partial differential equations-based model, it is difficult to get the data, so we do not validate our model with real data. But the model results are interesting and involve realistic assumptions.

Our work is divided in the following manner. In Section 2, we formulated our model, and using semigroup of operators, we show the existence and uniqueness of the solution. Section 3 is devoted to the stability analysis of disease-free equilibrium points. In Section 4, we study the optimal control problem. In Section 5, we solve our model numerically, and the last section is devoted to the concluding remarks.

2.   Model formulation

In this section, we formulate our model. We also show the existence and uniqueness of the solution to the model under consideration. We convert our problem into the abstract Cauchy problem to show the existence of the solution.

Let $S(a, t), E(a, t), I(a, t)$ and $R(a, t)$ denote the age and time-dependent population densities of individuals in the susceptible, exposed, infected, and recovered class. Whenever there is no confusion in the notation, we use the notation $S, E, I$ and $R$ for $S(a, t), E(a, t), I(a, t)$ and $R(a, t)$ respectively. We denote the age-dependent transmission coefficient by $c(a, b)$, which describes the interaction between susceptible and infected class, that is, $c(a_{1}, a_{2})S(a_{1}, t)I(a_{2}, t) da_{1} da_{2}$ is the number of susceptible individuals aged in $(a_{1}, a_{1}+da_{1})$ that contract the disease after interacting with infected individuals aged in $(a_{2}, a_{2}+da_{2})$. We assume that the functional form of the force of infection is given by

where $a_{m}$ is the maximum age which an individual can attain. Then the following system of partial differential equations (PDEs) describes the spread of disease

where $\varphi(a, t) = \int_{0}^{a_{m}}c(a, \sigma)I(\sigma, t) d \sigma$.

Let $X = \mathbb{R}^{4} \times L^{1}((0, a_{m}), \mathbb{R})\times L^{1}((0, a_{m}), \mathbb{R})\times L^{1}((0, a_{m}), \mathbb{R})\times L^{1}((0, a_{m}), \mathbb{R})$ and let us define the linear operator

with

Here it is clear that $\overline{D(A)}$ is not dense in $X$. Let us also define the nonlinear operator in the following manner

Our model in abstract form can be written as

where $v(t) = (0, 0, 0, 0, S(\cdot, t), E(\cdot, t), I(\cdot, t), V(\cdot, t))$ and $v_{0} = (0, 0, 0, 0, S_{0}, E_{0}, I_{0}, V_{0})$. Instead of solving (2.2), we consider the following integral equation

We have the following assumptions:

(i) $m, w_{0}, v, \sigma$ are positive functions and lies in $L_{+}^{\infty} ((0, a_{m}), \mathbb{R}).$

(ii) For $0 \le a \le a_{m}$, we assume that $0 < \sigma(a) < 1.$

Now, our task is to prove that $A$ is a Hille-Yosida operator and therefore it generates a $C_{0}$ semigroup on the closure of $D(A)$. Linearized system of (2.2) can be written as

where $v^{*}$ is steady state solution of (2.2) and

Here, $DF(v^{*})$ is a compact bounded linear operator on $X$. Let

where

Now, our task is to prove that $(A, D(A))$ is a Hille-Yosida operator. For

we have

This implies

Solving we get,

Now, taking $L^{1}$ norm we get

Therefore, we have

So, we conclude that $(A, D(A))$ is a Hille-Yosida operator.

Now, let us recall the following results.

Lemma 2.1. [Pazy ^[33], 1983] Let $X$ be a Banach space. Let operator $A$ with domain $D(A)$ be a Hille-Yosida operator on $X$ and $B$ be a bounded linear operator on $X$. Then $A+B$ is also a Hille-Yosida operator.

Lemma 2.2. [Pazy ^[33], 1983] Let $(A, D(A))$ be a Hille-Yosida operator, then its part $(A_{0}, D(A_{0}))$ generates a $C_{0}$ semigroup $\{ T_{0}(t)\; | \; t \ge 0 \}$ on $X_{0},$ where $X_{0} = (\overline{D(A)}, \| \cdot \|)$.

Using Lemmas 2.1 and 2.2, and proof given above, we conclude that the operator $A + DF(v^{*})$ is a Hille-Yosida operator and the part of $(A, D(A)$ and the part of $(A+DF(u^{*}))$ generate $C_{0}$ semigroups on $X_{0}$. So, the proof given above and Lemma 2.1, Lemma 2.2 prove the existence and uniqueness of mild solution to our model (for more details we refer to Pazy ^[33]). So, we have the following theorem:

Theorem 2.3. For any $v_{0} \in D_{+} = \{0 \} \times \{0 \} \times \{0 \} \times \{0 \} \times W^{1, 1}((0, a_{m}), \mathbb{R}_{+}) \times W^{1, 1}((0, a_{m}), \mathbb{R}_{+}) \times W^{1, 1}((0, a_{m}), \mathbb{R}_{+}) \times W^{1, 1}((0, a_{m}), \mathbb{R}_{+})$ system 2.4 has a unique continuous solution with initial data in $D_{+}$. Moreover, the function defined by $\Psi(t, v_{0}) = v(t, v_{0})$ is strongly continuous semigroup.

Remark 2.4. Since $\Psi(t, v_{0})$ is strongly continuous semigroup, there exists $\mathcal{M} \ge 1, \; \omega \ge 0$ such that $\| \Psi(t) \| \le \mathcal{M}e^{\omega t}$. So, if we assume $t \in [0, T]$ for some finite $T > 0$, then we have the boundedness of solution to the abstract Cauchy problem 2.2.

3.   Stability analysis

In this section, we discuss the stability analysis of disease-free steady state. Steady-state solutions play an important role to study the qualitative properties of the solution when the explicit form of the solution is not known. Steady-state equations for our model are given by

Disease free steady state $(S^{0}(a), 0, 0, 0, V^{0}(a))$ is given by

Now, our task is to derive expressions for $S^{0}$ which is independent of $V^{0}$ and similarly for $V^{0}$. Let

Now, Eqs (3.2) and (3.3) can be written as

Now, substituting value of $S^{0}$ from (3.5) into (3.4), we get

So, we have steady-state solutions in explicit form. From these solutions, it is clear that both $S^{0}$ and $V^{0}$ are non-negative. Let us define the threshold number in the following manner (for more details we refer to ^[27]):

The above defined threshold parameter may not be biologically relevant but plays an important role to study the stability of disease-free steady state. We can state the following result which relates the stability of disease-free steady state with the threshold parameter $R_{0}$.

Theorem 3.1. The disease free steady state is stable if $R_{0} < 1$ and unstable if $R_{0} > 1$.

Proof. For the local stability of the disease-free equilibrium point, firstly we will linearize our system around the disease-free steady state. Linearizing system (2.1) around the disease-free equilibrium, we get

We are looking for solutions of the form

Therefore, system (3.9) becomes

Solving these equations, we get

Substituting value of $\bar{E}(a)$ in $\bar{I}(a)$, we get

Observe that

In simplified form, it can be written as

Using (3.11), the characteristic equation can be written as

where

It is easy to see that $G(0) = R_{0}$. $G$ can also be written as

Taking derivative of $G$ with respect to $\lambda$, we get

Observe that $G$ is a decreasing function of $\lambda$ as $G'(\lambda) < 0$ and $\lim_{\lambda \to \infty}G(\lambda) = 0$. Let us assume that $\lambda = x_1 + ix_2$ is a root of equation $G(\lambda) = 1$. Then for $x_1 \ge 0$, we have

Thus the real part of eigenvalue $\lambda$ is negative if threshold parameter $R_0 < 1$. Therefore, the disease free equilibrium point is locally asymptotically stable if $R_0 < 1$ and unstable if $R_0 > 1$. □

4.   Optimal control problem

In this section, we discuss the optimal control strategy for our problem. Optimal control theory plays an important role in epidemiological models because it is very important to control the spread of disease in an optimal manner. Here, our control variable is vaccination effort and we want to minimize the cost of vaccination and also our aim is to reduce the number of infected individuals. We assume that $\varphi(a, t) = \psi(a, t)I(a, t)$. For this form of $\varphi(a, t)$, it is easy to construct operators $A$ and $F$ to show the existence and uniqueness of the solution. Also, we assume that $v$ depends on both age variable $a$ and time $t$. We consider the following optimal control problem:

$v(a, t)$ is the control variable and our aim is to minimize the number of infected individuals and the cost of implementing control.

Theorem 4.1. Let $(S^{v_{1}}, E^{v_{1}}, I^{v_{1}}, R^{v_{1}}, V^{v_{1}}), (S^{v_{2}}, E^{v_{2}}, I^{v_{2}}, R^{v_{2}}, V^{v_{2}})$ be solutions of system 4.1 corresponding to the control variables $v_{1}$ and $v_{2}$, respectively. Then for sufficiently small $T > 0$, the following estimates hold

where $\mathcal{C}_{01}$ is a positive constant which depends on the bound of $m, m_{1}, \gamma, \beta, \sigma, \psi, \omega_{0},$ and initial data. Moreover,

where $\mathcal{C}_{02}$ is a positive contant.

Proof. Using the method of characteristics, the system 4.1 can be solved as:

Let $\mathcal{Q} = (0, T) \times (0, a_{m})$, and assume that $(S^{v_{1}}, E^{v_{1}}, I^{v_{1}}, R^{v_{1}}, V^{v_{1}}), (S^{v_{2}}, E^{v_{2}}, I^{v_{2}}, R^{v_{2}}, V^{v_{2}})$ be solution of system (4.1) corresponding to the control variables $v_{1}$ and $v_{2}$ respectively. Then using the solution given by (4.2), we get

where $C_{11}, C_{22}, C_{33}$ are constants which depends on the bound of $m, \psi,$ and $\omega_{0}$. Similarly, we can consider the case when $a > t$, and we will obtain similar types of estimates. So, we can write

Using the same procedure for other state variables, we get

where $\mathcal{C}_{01}$ is a positive constant. Now, estimating the integral in the age variable only, we will get

where $\mathcal{C}_{02}$ is a positive constant that depends on model parameters and initial population distributions. More details about these types of estimates can be found in ^[31,32]. □

Let us define the following objective functional:

Here, $C_{1}, C_{2}$ and $C_{3}$ are weight factors and our aim is to minimize $J(v)$. We assume that the admissible control variable lies in the set $\mathcal{V}$ defined by

where we assume that $v \in L^{\infty}((0, a_{m}) \times (0, T))$ and $0 \le v \le v_{m}$. Now, our next result establish the existence of optimal control.

Theorem 4.2. Let $J: L^{1}(0, T) \to (-\infty, \infty]$ be defined by

Then to the control problem (4.1), there exists an optimal control $u^{*}$ such that

Proof. We have already shown the existence of unique solution to the system (4.1) which is uniformly bounded. Also control variable $v$ is uniformly bounded in $(0, a_{m})$. Therefore, there exist a sequence $v_{n} \in V$ such that

Because the function $C_{3}v^{2}$ is a lower semi-continuous function which is also convex, we have

Also, using the fact that the bound of $S, E, I, R$ and $V$ depends on the control variable $v$, we have $S_{n} = S(v_{n}), E_{n} = E(v_{n}), I_{n} = I(v_{n}), R_{n} = R(v_{n}), V_{n} = V(v_{n})$. Because $S_{n}, E_{n}, I_{n}, R_{n}$ and $V_{n}$ are uniformly bounded for each $n$. We have $S^{*} = S(v^{*}), E^{*} = E(v^{*}), I^{*} = I(v^{*}), R^{*} = R(v^{*})$ and $V^{*} = V(v^{*})$. Now,

From the relation $J(v^{*}) \le \inf_{v \in \mathcal{V}} J(v)$, we conclude that $v^{*}$ is optimal control which minimizes the objective functional. □

Now, our next task is to obtain the necessary conditions for our optimal control problem. Let us define $v^{\epsilon}(a, t) = v(a, t) + \epsilon g(a)$, where $0 < \epsilon < 1, \; g \in L^{\infty}(0, a_{m})$ and also called the variation function. Let us denote $S^{\epsilon} = S(v^{\epsilon}), E^{\epsilon} = E(v^{\epsilon}), I^{\epsilon} = I(v^{\epsilon}), R^{\epsilon} = R(v^{\epsilon})$ and $V^{\epsilon} = V(v^{\epsilon})$. Then we have the following system

Let us set the following notations for weak derivatives of $S, E, I, R$ and $V$:

$\bar{S}(a, t) = \lim_{\epsilon \to 0}\frac{S(v+\epsilon g)-S(v)}{\epsilon}, \bar{E}(a, t) = \lim_{\epsilon \to 0} \frac{E(v+\epsilon g)-E(v)}{\epsilon}, \bar{I}(a, t) = \lim_{\epsilon \to 0} \frac{I(v+\epsilon g)-I(v)}{\epsilon}, \bar{R}(a, t) = \lim_{\epsilon \to 0} \frac{R(v+\epsilon g)-R(v)}{\epsilon}, \bar{V}(a, t) = \lim_{\epsilon \to 0}\frac{V(v+\epsilon g)-V(v)}{\epsilon}$. Then we have the following system of differential equations

The directional derivative of $J$ with respect to $v$ in the direction of $g$ is given by

Now, our next task is to find the adjoint system. Let us define

Observe that

with the following conditions

Similarly we have

Now, we get the following adjoint system

Also, we have the following conditions

Using the same steps as we followed in 4.1, we can state the following theorem

Theorem 4.3. Let $(\lambda_{11}, \lambda_{21}, \lambda_{31}, \lambda_{41}, \lambda_{51})$ and $(\lambda_{12}, \lambda_{22}, \lambda_{32}, \lambda_{42}, \lambda_{52})$ be solution of adjoint system (4.19) with control variable $v_{1}$ and $v_{2}$ respectivey. Then

where $\mathcal{C}_{11}$ is a positive constant.

Now, we state and prove the following result, which gives the explicit form of optimal control.

Theorem 4.4. Let $(S^{*}, E^{*}, I^{*}, R^{*}, V^{*})$ and $(\lambda_{1}, \lambda_{2}, \lambda_{3}, \lambda_{4}, \lambda_{5})$ be state solutions and solutions to corresponding adjoint system respectively. Also, let $v^{*} \in \mathcal{V}$ be optimal control which minimizes (4.7). Then optimal control will be of the form

Proof. We know that

Now using the adjoint system, we have

Using (4.14, 4.15, 4.16, 4.17) and (4.18), we have

Because $J'(v) \ge 0$, so this implies in this case

By using the upper and lower bounds of control variable, we have

□

Remark 4.5. It is clear from the characterization of optimal control that it depends on the susceptible, infected, and vaccinated class. Although it depends on function $g$, we can always fix a particular choice of $g$.

Theorem 4.6. Let $M$ be a generic positive constant that depends on $\mathcal{C}_{02}, \mathcal{C}_{11}$ and the bound of vaccine vaning rate $\omega_{0}$, and transmission coefficient $\psi$. Then for sufficiently small $\frac{M}{2C_{3}}$, there exists a unique optimal controller $v^{*}$.

Proof. Let $(S^{v_{1}}, E^{v_{1}}, I^{v_{1}}, R^{v_{1}}, V^{v_{1}}), (S^{v_{2}}, E^{v_{2}}, I^{v_{2}}, R^{v_{2}}, V^{v_{2}})$ be solution of system (4.1) corresponding to the control variables $v_{1}$ and $v_{2}$ respectively, and $(\lambda_{11}, \lambda_{21}, \lambda_{31}, \lambda_{41}, \lambda_{51}), (\lambda_{12}, \lambda_{22}, \lambda_{32}, \lambda_{42}, \lambda_{52})$ be solution of adjoint system 4.19 with control variable $v_{1}$ and $v_{2}$ respectivey. Define the map $L: \mathcal{U} \to \mathcal{U}$ by

Then

where $M$ is a positive generic constant which depends on the bound of $\omega_{0}, g$ and $\psi$. Therefore, we have

where $M$ is a positive generic constant which also depends on $\mathcal{C}_{01}, \mathcal{C}_{02}$ and $T$. So, if $\frac{M}{2C_{3}} < 1$, then we have the existence of unique fixed point $v^{*}$. □

5.   Numerical simulation

In this section, we discuss numerical simulation results. We applied an explicit finite difference scheme to solve our model numerically. Since the explicit finite difference scheme is highly sensitive to the discretization size of independent variables, the discretization size is chosen in such a way that the scheme is convergent. A finite difference scheme is applied to solve the system of differential equations and for the adjoint system, we applied a backward finite difference scheme. {For the optimal control problem, we applied the forward-backward sweep method. The forward-backward sweep method for optimal control problems in PDEs combines the forward integration of the state equations with the backward integration of the adjoint equations. By iteratively adjusting the control variables based on the computed gradient, the method aims to find the optimal control that minimizes the cost functional while satisfying the given system of PDEs with initial and boundary data.} We use the following initial age distribution in our model:

Mortality rate is taken $\mu(a) = \frac{0.01}{a+0.0001}, \psi(a, t) = 0.02e^{-a}.$ Other parameters taken are $m_{1} = 0.2, q_{1} = 0.4, \gamma_{1} = 0.4, \gamma_{2} = 0.6, \gamma = 0.02, \beta = 1, \Lambda = 0.9$ and $g = 1$.

For the first two plots Figures 2(a), (b), we assume that a fixed proportion of individuals are vaccinated and observe the dynamics of infected individuals. So, for these figures we do not use the optimal vaccination strategy.

From Figure 2(a), (b), it is clear that increasing the vaccination rate will reduce the proportion of infected individuals. It is clear that increasing the vaccination waning rate will increase the population density of individuals who again become susceptible to disease. In Figure 3(a), (b), we have taken vaccine waning rates as 0 and 0.3 respectively. So, even if we are vaccinating a population, still there will be an increase in the population density of susceptible individuals due to imperfect vaccines.

Figure 4(a), (b) and Figure 5 show the dynamics of control variable for vaccine waning rate 0.1 and 0.01. As time progresses, the profile of the control variable attains a stable behavior. Initially, the control will be high and as time progresses, the control will take its minimum value.

In Section 3, we already prove that if threshold parameter $R_{0}$ is less than 1, then the steady state in which there are no individuals in the exposed or infected class is stable. In Figures 6 and 7, we have taken age between 10 and 50 and studied the stability of disease-free equilibrium point. We have chosen the parameters in such a way that the threshold number $R_{0} < 1$. These figures show that if threshold parameter $R_{0} < 1$, then the disease-free equilibrium point is stable. So, numerical simulations validate our theoretical results.

6.   Conclusions

Age-structured modeling is important when studying imperfect vaccination because it allows a more realistic representation of how infectious diseases spread in a population. Imperfect vaccination refers to situations where not everyone in a population is fully protected by a vaccine due to vaccine effectiveness or a decline in immunity over time.

In this work, we proposed a new age-structured SEIRV model, where we incorporate vaccine waning rate and imperfect vaccination in our model. The aim was to study the effect of leaky vaccines and to investigate the optimal approach to contain the disease spread. The Well-posedness of the proposed model is shown using the semigroup of operators. We convert our model into an abstract Cauchy problem and then show that the operator in a homogeneous problem is the Hille-Yosida operator. Finally, the basic reproduction number is defined as a threshold parameter to study the stability analysis of the disease-free equilibrium point.

Vaccination strategy and optimal control problems with vaccination effort as a control variable are also discussed. The optimality conditions are derived for the optimal control problem. With the help of the adjoint system, optimal control is characterized. Numerical simulations are also added to establish theoretical results. An explicit finite difference scheme is used to solve the model numerically. The discretization size is chosen so that the scheme is convergent. The stability of disease-free steady state is shown using numerical simulation. Using the forward-backward sweep method, we plot the optimal control variable with respect to age and time variables.

Age-structured compartmental models are a valuable tool for understanding and analyzing the effects of imperfect vaccination. Although getting real data for verification is difficult, the model assumptions and outcomes are quite realistic. This study can help the public health policymakers optimally vaccinate a population and reduce the impact of the disease with a low cost of vaccination. The threshold number $R_0$ involves some parameters, and by increasing or decreasing these parameters, accordingly, we can eradicate the infectious disease. By choosing different vaccine waning rates, the dynamics of individuals in different compartments can be studied using the numerical simulation results. We also encourage other researchers to explore this area. {Moreover, the age-structured SEIR model with stochastic perturbation can be investigated. More general size-structured variables can also be considered for epidemiological models. Finally, It is also natural to consider the delay in age and size-structured population models as it represents the maturation period. Therefore, age-structured compartmental models with delay can also be explored.

Use of AI tools declaration

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

Acknowledgement

The authors would like to thank the anonymous reviewers for their valuable comments and suggestions that helped improve this manuscript's quality.

A. Tridane is supported by UAEU UPAR grant number 12S125.

Conflict of interest

The authors declare there is no conflict of interest.