Studying the Impact of Vaccination Strategy and Key Parameters on Infectious Disease Models

Tahmineh Azizi, Bacim Alali

Department of Mathematics, Kansas State University, Manhattan, KS, USA.

DOI: 10.4236/ojop.2020.93007                

Abstract

In the current work, we study two infectious disease models and we use nonlinear optimization and optimal control theory which helps to find strategies towards transmission control and to forecast the international spread of the infectious diseases. The relationship between epidemiology, mathematical modeling and computational tools lets us to build and test theories on the development and fighting with a disease. This study is motivated by the study of epidemiological models applied to infectious diseases in an optimal control perspective. We use the numerical methods to display the solutions of the optimal control problems to find the effect of vaccination on these models. Finally, global sensitivity analysis LHS Monte Carlo method using Partial Rank Correlation Coefficient (PRCC) has been performed to investigate the key parameters in model equations. This present work will advance the understanding about the spread of infectious diseases and lead to novel conceptual understanding for spread of them.

Keywords

Optimal Control, S-I-R Model, S-E-I-R Model, LHS Monte Carlo Method, Fourth Order Runge-Kutta

1. Introduction

Recently, due to the fast spread of pandemic diseases, mathematical modeling in the field of epidemiology has attracted many scientists in different areas. Many mathematical models have been developed to describe the transmission of communicable diseases [1] [2] [3]. These mathematical models describe the mechanisms of infectious diseases and also, they are helpful to analyze the effect of public health interventions to control the spreading of diseases. In mathematical perspective, we describe biological systems by converting them into mathematical and theoretical framework with biological parameters and then using computer code to solve the model system computationally to predict the future of infectious diseases, one needs to study the behavior of each individual which plays a key role to understand the behavior epidemiology of infectious diseases [4].

One of the new approaches in modeling dynamic systems is the theory of optimal control. For the first time, R. E. Bellman introduced a new method to solve dynamic systems by using the principle of optimality which reduces significantly the computation of the optimal controls [5]. In optimal control (OC) theory, for a dynamic system, we define a control problem and its state trajectories over a period of time to minimize a performance index [6]. In optimal control theory, the problem of determining the control would be turned to an extension of the calculus of variations [7]. One of the most interesting applications of the calculus of variations was in the Hamilton’s principle or the Principle of Least Action. The Russian mathematician Lev S. Pontryagin and his colleagues V. G. Boltyanskii, R. V. Gamkrelidz and E. F. Misshchenko generalized the calculus of variations to optimal control theory by proposing the Pontryagin Maximum Principle [8] which defines appropriate conditions for optimization problems with differential equations as constraints. OC can be used for the problems where the calculus of variations is not applicable, such as the problems which include constraints on the derivatives of functions [9]. With increasing the number of variables and parameters of system, optimal control problems cannot be solved analytically and one may need to apply numerical methods.

To model a dynamic system, we usually use a set of ordinary differential equations. A system of ODEs for $t_0\le t\le t_f$ can be described by

$\stackrel{\dot{}}{y}=\left[\begin{array}{c}{\stackrel{\dot{}}{y}}_1\\ {\stackrel{\dot{}}{y}}_2\\ ⋮\\ {\stackrel{\dot{}}{y}}_n\end{array}\right]=\left[\begin{array}{c}{f}_1\left(y_1\left(t\right),\cdots ,y_n\left(t\right),t\right)\\ f_2\left(y_1\left(t\right),\cdots ,y_n\left(t\right),t\right)\\ ⋮\\ f_n\left(y_1\left(t\right),\cdots ,y_n\left(t\right),t\right)\end{array}\right]$

Based on how the conditions at the endpoints of the domain are specified, we classify an ODE solving problem into initial value problems (IVP) and boundary value problems (BVP). For an initial-value problem, all the conditions are specified at the initial point. For a boundary-value problem, the conditions are needed for both initial and final points. There are many numerical methods to solve initial value problems, such as Euler, Runge-Kutta or adaptive methods and boundary value problems, such as shooting methods [10] [11].

Euler method is the most common used single-step method. In this discretization technique, for differential equation $\stackrel{\dot{}}{x}=f\left(x\left(t\right)\mathrm{,}t\right)$, we can make a convenient approximation of this:

$x_{n+1}\simeq x_n+hf\left(x\left(t_n\right)\mathrm{,}{t}_n\right):$

The approximation $x_{n+1}$ of $x\left(t\right)$ at the point $t_{n+1}$ has an error of order h². There exists a trade-off between accuracy and complexity of calculation which depends heavily on the chosen value for h. As h is decreasing, the calculation would be longer however more exact. One of the disadvantages of this method is for many higher order systems. It is very difficult to have an effective Euler approximation. Thus, we need to use more accurate and elaborate methods and one of them is the Runge-Kutta method.

Runge-Kutta method is a multiple-step method. In this technique, we obtain the solution at time $t_{k+1}$ from the values $t_{j-k}\mathrm{,}\cdots \mathrm{,}{t}_k$ and j is the number of steps. To approximate a differential equation of the form $\stackrel{\dot{}}{x}=f\left(x\left(t\right)\mathrm{,}t\right)$, we can use the second order Runge-Kutta method

$x_{n+1}\simeq x_n+\frac{h}{2}\left[f\left(x_n\left(t\right)\mathrm{,}{t}_n\right)+f\left(x_{n+1}\left(t\right)\mathrm{,}{t}_{n+1}\right)\right]\mathrm{<mo>\; ;\; </mo>}$

or the fourth order Runge-Kutta method

$x_{n+1}\simeq x_n+\frac{h}{6}\left(k_1+2k_2+2k_3+k_4\right)$

where

$k_1=f\left(x\left(t\right)\mathrm{,}t\right)$

$k_2=f\left(x\left(t\right)+\frac{h}{2}{k}_1\mathrm{,}t+\frac{h}{2}\right)$

$k_3=f\left(x\left(t\right)+\frac{h}{2}{k}_2\mathrm{,}t+\frac{h}{2}\right)$

$k_4=f\left(x\left(t\right)+h{k}_3\mathrm{,}t+h\right)$

For the second and fourth order Runge-Kutta method, the approximation $x_{n+1}$ of $x\left(t\right)$ at the point $t_{n+1}$ has an error of order $h_3$ and $h_5$.

In this research, we study the most basic epidemiological models S-I-R model (composed by Susceptible-Infected-Recovered) and S-E-I-R model (Susceptible-Exposed-Infected-Recovered). For these models, we develop some analytical results that are useful in understanding of simple epidemic diseases. We continue this study by proposing the equivalent optimal control problems of the mentioned epidemic models and we numerically solve them using the backward-forward sweep method with fourth order Runge-Kutta. Finally, we perform global sensitivity analysis by LHS Monte Carlo method using PRCC to identify the key parameters that contribute most significantly to the spread or control of the infectious diseases.

2. Kermack-McKendrick SIR Epidemic Model (S-I-R Model)

Recently, due to the fast spread of pandemic diseases, mathematical modeling in the field of epidemiology has attracted many scientists in different areas. Many mathematical models have been developed to describe the transmission of communicable diseases and among these models, the classical Kermack-McKendrick SIR epidemic model builds the basic skeleton of all of them [12].

S := Susceptible (People who could potentially catch the disease)

I := Infective (People who currently have the disease)

R := Removed (People recovered or have died)

Assumptions:

1) Total population remains constant

2) Rate of increase in the infectives is proportional to the contact between susceptible and infective

3) Removal rate (death rate is constant)

Using these assumptions, the classical S-I-R model has the following form:

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=-\beta IS+\delta R\\ \frac{\text{d}I}{\text{d}t}=\beta IS-\gamma I\\ \frac{\text{d}R}{\text{d}t}=\gamma I-\delta R\end{array}$ (2.1)

where, $\beta$ demonstrates rate of infection, $\gamma$ implies to rate of recovery and $\delta$ represents rate of immunity loss. If $\delta =0$, we assume a model without immunity loss. In the first equation of system (2.1), susceptible S decreases according to the number of contacts between infective I and susceptible S. Therefore, because of decreasing the rate of change of susceptible over time, in the first equation we get $-\beta IS$. The rate of change of infective I increases by IS and decreases by $\gamma I$. The term $\beta IS$ has been added to the second equation of system (2.1) which is due to the increasing the contact between S and I. The negativity of $\gamma I$ is showing decreasing the rate of change in infective I by moving to the next stage which is recovered or died. The term $\gamma I$ has been added t the third equation which means that the rate of changing the recovered R is increasing by this factor. The time-evolution of system (2.1) over 300 days have been demonstrated in Figures 1-4.

Will the Disease Spread? What Is the Max Number of Infectives I_max? How Many People Catch the Disease?

To answer these questions consider the following general S-I-R model:

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=-\beta IS\\ \frac{\text{d}I}{\text{d}t}=\beta IS-\gamma I\\ \frac{\text{d}R}{\text{d}t}=\gamma I\end{array}$ (2.2)

At the start of outbreak we have $S=S_0$, $I=I_0$ and $R=0$. Total population size remains constant during epidemic; therefore, the rate of change of $S+I+R$ must be zero:

$\frac{\text{d}}{\text{d}t}\left(S+I+R\right)=0,S+I+R=S_0+I_0$ (2.3)

Figure 1. The SIR schematic model for system (2.1). S := Susceptible Compartment, I := Infective Compartment, R := Removed Compartment.

Figure 2. The time-evolution of disease over 300 days $\beta =5\times {10}^{-9}$, $\gamma =0.12$, $\delta =0.016$.

Figure 3. The time-evolution of disease over 300 days $\beta =5\times {10}^{-9}$, $\gamma =0.12$, $\delta =0.0$.

To find out if the disease will spread, we need to check that

$\frac{\text{d}I}{\text{d}t}=I\left(\beta S-\gamma \right)<0,S\le S_0$

Therefore, if $S_0>\frac{\gamma }{\beta }=\frac{1}{q}$, then disease will spread. Here, $\frac{1}{q}$ is the contact ration

Figure 4. The time-evolution of disease over 300 days $\beta =5\times {10}^{-9}$, $\gamma =0.07$, $\delta =0.0$.

which is the fraction of population that comes to contact with individual during the period of infectious. However, if the reproductive number or the ratio number

$R_0=\frac{\beta S_0}{\gamma }>1$, we have epidemic. This ratio represents the number of secondary

infection in the population caused by initial primary infection, i.e. how many other people get the disease.

To find the maximum number of infectives or $I_{\max }$, we combine $\frac{\text{d}S}{\text{d}t}$ and $\frac{\text{d}I}{\text{d}t}$ :

$\frac{\text{d}I}{\text{d}S}=\frac{\beta IS-\gamma I}{-\beta IS}=-1+\frac{\gamma }{\beta S}=-1+\frac{1}{qS}$

Assuming

$I+S-\frac{1}{q}\ln S=I_0+S_0-\frac{1}{q}\ln S_0$ (2.4)

Then

$I_{\max }=I_0+S_0-\frac{1}{q}\left(1+\ln \left(q{S}_0\right)\right)$

Here, $I_{\max }$ represents the maximum number of people who have the disease at a given time. For COVID-19, or similar worldwide diseases the value for q (contact parameter) is high since disease easily transmits. When q is large, it means that the number of people get infected is a lot.

To reduce the reproduction rate, one can reduce the number of susceptible, $S_0$. One way to decrease the number of susceptible is using vaccination which is a common method to eradicate of infectious diseases. Vaccination can go further than being used for just individuals, but it can be beneficial in large scale communities by preserving the effective reproduction rate below the level which would allow an epidemic to spread. However, an epidemic can start and spread very quickly if the reproduction rate rises beyond the critical value for an epidemic [13].

To find out how many people catch the disease, based on the first assumption, the total population is constant and to end the disease, the number of infected need to go down to zero (end of out break):

$S+I+R=S_0+I_0$

and

$R\left(end\right)=-S\left(end\right)+I_0+S_0$

Here, $S\left(end\right)$ is unknown. From (2.4), we have

$S\left(end\right)-\frac{1}{q}\ln \left(S\left(end\right)\right)=I_0+S_0-\frac{1}{q}\ln S_0$

The graph of $S\left(end\right)$ is decreasing and shows at small value of $S\left(end\right)$ and larger q, we have larger value for $R\left(end\right)$.

3. The S-E-I-R Model

An SIR model is an epidemiological model that represents the number of people infected with a contagious illness in a closed population over time. In another word, there are some other important infections which include a significant latency or incubation period during which individuals have been infected but are not yet infectious themselves (for example this latency period is zero for cold). During incubation period the individual is exposed. See Figure 5.

Here, we write the total population as $N=S+E+I+R$. So, the S-E-I-R model has the form

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=\Lambda -\beta IS-\delta S\\ \frac{\text{d}E}{\text{d}t}=\beta IS-\epsilon E-\delta E\\ \frac{\text{d}I}{\text{d}t}=\epsilon E-\gamma I-\delta I\\ \frac{\text{d}R}{\text{d}t}=\gamma I-\delta R\end{array}$ (3.1)

Figure 5. The transport diagram for S-E-I-R model (2.1). S := Susceptible Compartment, E := Exposed Compartment, I := Infective Compartment, R := Removed Compartment.

where

S := Susceptible (People who could potentially catch the disease)

E := Exposed (People who are infected but are not yet infectious)

I := Infective (People who currently have the disease)

R := Removed (People recovered or have died)

$\delta$ := Constant death rate

$\Lambda =\mu \times N$ := Constant influx of new susceptible ( $\mu$ Constant birth rate)

$\epsilon$ := Latency transfer rate to infectious

$\gamma$ := Recovery rate of infectious

$\beta IS$ := The bilinear (mass action) incidence.

For simplicity, we assume that the death rates are equal ${\delta }_S={\delta }_E={\delta }_I={\delta }_R$.

If we have $S>0$ and $E=I=R=0$, we have a disease free population or disease free equilibrium, which means that there is no disease. To find disease or endemic equilibrium point, we look for a feasible region $\Sigma$ such that:

$\frac{\text{d}N}{\text{d}t}>0\to \frac{\text{d}\left(S+E+I+R\right)}{\text{d}t}=\frac{\text{d}S}{\text{d}t}+\frac{\text{d}E}{\text{d}t}+\frac{\text{d}I}{\text{d}t}+\frac{\text{d}R}{\text{d}t}>0$

Therefore, from (5) we have

$\begin{array}{l}\left(\Lambda -\beta IS-\delta S\right)+\left(\beta IS-\epsilon E-\delta E\right)+\left(\epsilon E-\gamma I-\delta I\right)+\left(\gamma I-\delta R\right)\\ =\Lambda -\delta S-\delta E-\delta I-\delta R=\Lambda -\delta \left(S+E+I+R\right)=\Lambda -\delta N>0\end{array}$

Thus

$\frac{\text{d}N}{\text{d}t}=\Lambda -\delta N\ge 0\Rightarrow \Lambda \ge \delta N\iff \frac{\Lambda }{\delta }\ge N=S+E+I+R$

Therefore, the feasible region $\Sigma$ would be:

$\Sigma =\left\{\left(S\mathrm{,}E\mathrm{,}I\mathrm{,}R\right)\in {ℝ}^4\mathrm{<mo>\; |\; </mo>}S+E+I+R\le \frac{\Lambda }{\delta }\right\}$

From equation three we have,

$\epsilon =\left(\delta +\gamma \right)I\Rightarrow I^{*}=\left(\frac{\epsilon }{\delta +\gamma }\right)E$

From equation four,

$\gamma I-\delta R=0\Rightarrow R^{*}=\frac{\gamma }{\delta }\left(\frac{\epsilon }{\delta +\gamma }\right)E$

1) Case 1: If $E^{*}=0$ No Exposed. So, $I^{*}=0$ and $R^{*}=0$.

From equation one,

$0=\Lambda -\delta S^{*}\Rightarrow S^{*}=\frac{\Lambda }{\delta }$

Therefore, the diseases free equilibrium would be:

$P_0=\left(\frac{\Lambda }{\delta },0,0,0\right)$

2) Case 2: If $E^{\mathrm{<mo>\; *\; </mo>}}\ne 0$, then $I^{\mathrm{<mo>\; *\; </mo>}}\ne 0$ and $R^{\mathrm{<mo>\; *\; </mo>}}\ne 0$.

If we add equation two and three, we get

$\left(\beta I^{*}{S}^{*}-\epsilon E^{*}-\delta E^{*}\right)+\left(\epsilon E^{*}-\gamma I^{*}-\delta I^{*}\right)=0+0=0$

$\Rightarrow \beta I^{\mathrm{<mo>\; *\; </mo>}}{S}^{\mathrm{<mo>\; *\; </mo>}}-\left(\delta +\gamma \right)I^{\mathrm{<mo>\; *\; </mo>}}-\delta E^{\mathrm{<mo>\; *\; </mo>}}=0$

$\Rightarrow \beta I^{\mathrm{<mo>\; *\; </mo>}}{S}^{\mathrm{<mo>\; *\; </mo>}}-\left(\epsilon +\delta \right)E^{\mathrm{<mo>\; *\; </mo>}}=0$

$\Rightarrow \beta \left(\frac{\epsilon }{\delta +\gamma }\right)E^{\mathrm{<mo>\; *\; </mo>}}{S}^{\mathrm{<mo>\; *\; </mo>}}-\left(\epsilon +\delta \right)E^{\mathrm{<mo>\; *\; </mo>}}=0$

$\Rightarrow \left[\beta \left(\frac{\epsilon }{\delta +\gamma }\right)S^{\mathrm{<mo>\; *\; </mo>}}-\left(\epsilon +\delta \right)\right]E^{\mathrm{<mo>\; *\; </mo>}}=0$

$\Rightarrow \left(\frac{\epsilon \beta }{\delta +\gamma }\right)S^{\mathrm{<mo>\; *\; </mo>}}-\left(\epsilon +\delta \right)=0$

$\Rightarrow S^{\mathrm{<mo>\; *\; </mo>}}=\frac{\left(\epsilon +\delta \right)\left(\delta +\gamma \right)}{\epsilon \beta }$

General replication number $R_0$ is the number of new cases any single infected individual is going to create and produce or infect susceptible. To find $R_0$, at

equilibrium we have $S^{\mathrm{<mo>\; *\; </mo>}}\ge \frac{\Lambda }{\delta }$. For $S^{\mathrm{<mo>\; *\; </mo>}}=\frac{\left(\epsilon +\delta \right)\left(\delta +\gamma \right)}{\epsilon \beta }$, we have:

$\frac{\left(\epsilon +\delta \right)\left(\delta +\gamma \right)}{\epsilon \beta }\ge \frac{\Lambda }{\delta }\iff 1\ge \frac{\Lambda \epsilon \beta }{\delta \left(\epsilon +\delta \right)\left(\delta +\gamma \right)}=R_0$

1) $R_0\le 1$ $\Rightarrow$ Disease Free equilibrium $P_0=\left(\frac{\Lambda }{\delta },0,0,0\right)$ ; we can control and there is no disease.

2) $R_0>1$ $\Rightarrow$ Endemic equilibrium $P^{*}=\left(S^{*},E^{*},I^{*},R^{*}\right)>0$.

The time-evolution of system (5) over 300 days have been demonstrated in Figure 6.

Figure 6. The time-evolution of system (3.1) over 300 days $\beta =5\times {10}^{-9}$, $\gamma =0.07$, $\delta =1/60$ and $\mu =1/50$.

4. Optimal Control Problem

A general optimal control (OC) problem needs a cost functional $\left(J\left[x\left(t\right)\mathrm{,}u\left(t\right)\right]\right)$, a set of state variables $\left(x\left(t\right)\in X\right)$, a set of control variables $\left(u\left(t\right)\in U\right)$ in a time t, with $t_0\le t\le t_f$. The main goal is finding a piecewise continuous control $u\left(t\right)$ and the associated state variable $x\left(t\right)$ to maximize a given objective functional.

Definition 4.1 (Basic OC Problem in Lagrange formulation). An OC problem is in the form

$\underset{u}{\max }=J\left[x\left(t\right),u\left(t\right)\right]={\displaystyle {\int }_{t_0}^{t_f}}f\left(t,x\left(t\right),u\left(t\right)\right)\text{d}t$

$\text{s}\text{.t}\stackrel{\dot{}}{x}\left(t\right)=g\left(t\mathrm{,}x\left(t\right)\mathrm{,}u(\; t\; )\right)$

$x\left(t_0\right)=x_0$

$x\left(t_f\right)$ could be free, which means that the value of $x\left(t_f\right)$ is unrestricted, or could be fixed, i.e., $x\left(t_f\right)=x$ [14].

We consider f and g to be continuously differentiable functions. We suppose that the control set U is a Lebesgue measurable function. Therefore, as long as the controls will always be piecewise continuous, the associated states will be piecewise differentiable.

We can change the maximization problem to a minimization problem by making the cost functional negative:

$\min \left\{J\right\}=-\max \left\{J\right\}$

Definition 4.2 (Bolza formulation). The Bolza formulation of the OC problem can be defined as

$\underset{u}{\max }=J\left[x\left(t\right),u\left(t\right)\right]=\Phi \left(t_0,x\left(t_0\right),t_f,x\left(t_f\right)\right)+{\displaystyle {\int }_{t_0}^{t_f}}f\left(t,x\left(t\right),u\left(t\right)\right)\text{d}t$

$\text{s}\text{.t}\stackrel{\dot{}}{x}\left(t\right)=g\left(t\mathrm{,}x\left(t\right)\mathrm{,}u(\; t\; )\right)$

$x\left(t_0\right)=x_0$

where $\Phi$ is a continuously differentiable function [15].

Definition 4.3 (Mayer formulation). [16] The Mayer formulation of the OC problem can be defined as

$\underset{u}{\max }=J\left[x\left(t\right),u\left(t\right)\right]=\Phi \left(t_0,x\left(t_0\right),t_f,x\left(t_f\right)\right)$

$\text{s}\text{.t}\stackrel{\dot{}}{x}\left(t\right)=g\left(t\mathrm{,}x\left(t\right)\mathrm{,}u(\; t\; )\right)$

$x\left(t_0\right)=x_0$

4.1. Pontryagin’s Maximum Principle

Pontryagin proposed the idea of adjoint functions to append the differential equation to the objective functional which was one of the most important results of Mathematics in the 20th century and illustrates the necessary conditions to find the optimal control. Similar to Lagrange multipliers in multivariate calculus, Adjoint functions append constraints to the function of several variables to be maximized or minimized [7].

Definition 4.4 (Hamiltonian). Consider the OC problem in definition (4.1). The function

$H\left(t,x\left(t\right),u\left(t\right),\lambda \left(t\right)\right)=f\left(t,x\left(t\right),u\left(t\right)\right)+\lambda \left(t\right)g\left(t,x\left(t\right),u(\; t\; )\right)$

is called Hamiltonian function and $\lambda$ is the adjoint variable.

Theorem 4.5 (Pontryagin’s Maximum Principle). [8] [17] Let $u^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)$ and $x^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)$ be optimal for problem in definition (4.1), then there exists a piecewise differentiable adjoint variable $\lambda \left(t\right)$ such that

$H\left(t\mathrm{,}{x}^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}u\left(t\right)\mathrm{,}\lambda \left(t\right)\right)\le H\left(t\mathrm{,}{x}^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{u}^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}\lambda (\; t\; )\right)$

for all controls u at each time t, where H is the Hamiltonian previously defined and

${\lambda }^{\prime }\left(t\right)=\frac{\partial H\left(t\mathrm{,}{x}^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{u}^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}\lambda \left(t\right)\right)}{\partial x}$

$\lambda \left(t_f\right)=0$

The last condition, $\lambda \left(t_f\right)=0$, called transversality condition, is only used when the OC problem does not have terminal value in the state variable, i.e., $x\left(t_f\right)$ is free.

This Pontryagin’s Maximum Principle converts the problem of finding a control which maximizes the objective functional subject to the state ODE and initial condition into the problem of optimizing the Hamiltonian pointwise. Therefore, with this adjoint equation and Hamiltonian, we have

$\frac{\partial H}{\partial u}=0$

at $u^{\mathrm{<mo>\; *\; </mo>}}$ for each t, meaning that the Hamiltonian has a critical point and we call this condition as optimality condition. Therefore, to find the necessary conditions, we do not need to calculate the integral in the objective functional and we only use the Hamiltonian.

4.2. Existence of a Finite Objective Functional Value at the Optimal Control and State Variables

Theorem 4.6. [18] [19] Consider

$\underset{u}{\max }=J\left[x\left(t\right),u\left(t\right)\right]={\displaystyle {\int }_{t_0}^{t_f}}f\left(t,x\left(t\right),u\left(t\right)\right)\text{d}t$

$\text{s}\text{.t}\stackrel{\dot{}}{x}\left(t\right)=g\left(t,x\left(t\right),u(\; t\; )\right)$

$x\left(t_0\right)=x_0$

Suppose that $f\left(t\mathrm{,}x\mathrm{,}u\right)$ and $g\left(t\mathrm{,}x\mathrm{,}u\right)$ are both continuously differentiable functions in their three arguments and concave in x and u. Suppose $u^{\mathrm{<mo>\; *\; </mo>}}$ is a control with associated state $x^{\mathrm{<mo>\; *\; </mo>}}$, and $\lambda$ a piecewise differentiable function, such that $u^{\mathrm{<mo>\; *\; </mo>}}$, $x^{\mathrm{<mo>\; *\; </mo>}}$ and $\lambda$ together satisfy on $t_0\le t\le t_f$ :

$f_u+\lambda g_u=0,$

${\lambda }^{\prime }=-\left(f_x+\lambda g_x\right)\mathrm{,}$

$\lambda \left(t_f\right)=\mathrm{0,}$

$\lambda \left(t\right)\ge 0.$

Then for all controls u, we have $J\left(u^{\mathrm{<mo>\; *\; </mo>}}\right)\ge J(\; u\; )$

Based on how the conditions at the endpoints of the domain are specified, we classify an ODE solving problem into initial value problems (IVP) and boundary value problems (BVP). For an initial-value problem, all the conditions are specified at the initial point. For a boundary-value problem the conditions are needed for both initial and final points. There are many numerical methods to solve initial value problems such as Euler, Runge-Kutta or adaptive methods and boundary value problems, such as shooting methods [10] [11].

Numerical methods for solving OC problems started from the 1950s with the works of Bellman [10]. We can divide these methods into two major groups: direct methods and indirect methods. Indirect methods indirectly solve the problem by converting the optimal control problem to a boundary-value problem, using the PMP. However, direct method solves the OC by transcribing an infinite-dimensional optimization problem to a finite-dimensional optimization problem.

5. OC Problem for S-I-R Model

In this section, we present an optimal control (OC) problem to study the dynamics of S-I-R model, using a vaccination process (u) as a measure to control the disease. Let $x_1$ represents the susceptible population, $x_2$ represents the proportion of population that is infected and $x_3$ represents the proportion of population that is recovered or dead. The optimal control problem can be defined as:

$\underset{u}{\min }=J\left[x\left(t\right),u\left(t\right)\right]={\displaystyle {\int }_{t_0}^{t_f}}\left(x_2+u^2\right)\text{d}t$ (5.1)

$\text{s}\text{.t}\frac{\text{d}{x}_1}{\text{d}t}=-\beta x_1{x}_2+\delta x_3-u{x}_1$ (5.2)

$\frac{\text{d}{x}_2}{\text{d}t}=\beta x_1{x}_2-\gamma x_2$ (5.3)

$\frac{\text{d}{x}_3}{\text{d}t}=\gamma x_2-\delta x_3$ (5.4)

$x\left(t_0\right)=\left(x_1\left(0\right),x_2\left(0\right),x_3\left(0\right)\right)$ (5.5)

with $x\left(t\right)=\left(x_1\left(t\right),x_2\left(t\right),x_3\left(t\right)\right)$ and $\lambda \left(t\right)=\left({\lambda }_1\left(t\right),{\lambda }_2\left(t\right),{\lambda }_3\left(t\right)\right)$, with initial conditions $x_1\left(0\right)=6\times {10}^7$, $x_2\left(0\right)={10}^7$, $x_3\left(0\right)=10$ and the parameters $\beta =5\times {10}^{-9}$, $\gamma =0.12$, $\delta =1/60$.

Let consider the problem (5.1) and constraints (5.2)-(5.4). With $x\left(t\right)=\left(x_1\left(t\right),x_2\left(t\right),x_3\left(t\right)\right)$ and $\lambda \left(t\right)=\left({\lambda }_1\left(t\right),{\lambda }_2\left(t\right),{\lambda }_3\left(t\right)\right)$, the Hamiltonian of this problem can be written as

$\begin{array}{c}H\left(t,x\left(t\right),u\left(t\right),\lambda \left(t\right)\right)=A{x}_2+u^2+{\lambda }_1\left(-\beta x_1{x}_2+\delta x_3-u{x}_1\right)\\ +{\lambda }_2\left(\beta x_1{x}_2-\gamma x_2\right)+{\lambda }_3\left(\gamma x_2-\delta x_3\right)\end{array}$

A is a weight parameter describing the comparative importance of the two terms in the functional. Using the PMP the optimal control problem can be studied with the state variables

${\stackrel{\dot{}}{x}}_1=-\beta x_1{x}_2+\delta x_3-u{x}_1$

${\stackrel{\dot{}}{x}}_2=\beta x_1{x}_2-\gamma x_2$

${\stackrel{\dot{}}{x}}_3=\gamma x_2-\delta x_3$

The adjoint variables are:

${\stackrel{\dot{}}{\lambda }}_1={\lambda }_1\left(u+\beta x_2\right)+{\lambda }_2\beta x_2$

${\stackrel{\dot{}}{\lambda }}_2=-A+{\lambda }_1\beta x_1-{\lambda }_2\left(\beta x_1-\gamma \right)-{\lambda }_3\gamma$

${\stackrel{\dot{}}{\lambda }}_3={\lambda }_3\delta -{\lambda }_1\delta$

with transversality conditions $\lambda \left(t_f\right)=0$. Figure 7 demonstrates the optimal curves for the states variables and optimal control corresponding to S-I-R model (2.1).

6. OC problem for S-E-I-R Model

In this section, we present an optimal control (OC) problem to study the dynamics of S-E-I-R model, using a vaccination process (u) as a measure to control the disease. Let $x_1$ represents the susceptible population, $x_2$ represents the

Figure 7. Solutions of optimal control problem for S-E-I-R model (2.1). u := Vaccination related variable, S := Susceptible Population, I := Infective Population, R := Removed Population.

proportion of population that is in the incubation period, $x_3$ represents the proportion of population that is infected and $x_4$ represents the proportion of population that is recovered or dead. The optimal control problem can be defined as:

$\underset{u}{\min }=J\left[x\left(t\right),u\left(t\right)\right]={\displaystyle {\int }_{t_0}^{t_f}}\left(x_3+u^2\right)\text{d}t$ (6.1)

$\text{s}\text{.t}\frac{\text{d}{x}_1}{\text{d}t}=\Lambda -\beta x_1{x}_3-\delta x_1-u{x}_1$ (6.2)

$\frac{\text{d}{x}_2}{\text{d}t}=\beta x_1{x}_3-\epsilon x_2-\delta x_2$ (6.3)

$\frac{\text{d}{x}_3}{\text{d}t}=\epsilon x_2-\gamma x_3-\delta x_3$ (6.4)

$\frac{\text{d}{x}_4}{\text{d}t}=\gamma x_3-\delta x_4$ (6.5)

$x\left(t_0\right)=\left(x_1\left(0\right),x_2\left(0\right),x_3\left(0\right),x_4\left(0\right)\right)$ (6.6)

with initial conditions $x_1\left(0\right)=6\times {10}^7$, $x_2\left(0\right)={10}^7$, $x_3\left(0\right)=10$, $x_4\left(0\right)=1$ and the parameters $\beta =5\times {10}^{-9}$, $\gamma =0.12$, $\delta =1/60$.

Let consider the problem (6.1) and constraints (6.2)-(6.5). With $x\left(t\right)=\left(x_1\left(t\right)\mathrm{,}{x}_2\left(t\right)\mathrm{,}{x}_3\left(t\right)\mathrm{,}{x}_4\left(t\right)\right)$ and $\lambda \left(t\right)=\left({\lambda }_1\left(t\right)\mathrm{,}{\lambda }_2\left(t\right)\mathrm{,}{\lambda }_3\left(t\right)\mathrm{,}{\lambda }_4\left(t\right)\right)$, the Hamiltonian of this problem can be written as

$\begin{array}{c}H\left(t\mathrm{,}x\left(t\right)\mathrm{,}u\left(t\right)\mathrm{,}\lambda \left(t\right)\right)=A{x}_3+u^2+{\lambda }_1\left(\Lambda -\beta x_1{x}_3-\delta x_1-u{x}_1\right)\\ +{\lambda }_2\left(\beta x_1{x}_3-\epsilon x_2-\delta x_2\right)+{\lambda }_3\left(\epsilon x_2-\gamma x_3-\delta x_3\right)\\ +{\lambda }_4\left(\gamma x_3-\delta x_4\right)\end{array}$

A is a weight parameter describing the comparative importance of the two terms in the functional. Using the PMP the optimal control problem can be studied with the state variables

${\stackrel{\dot{}}{x}}_1=\Lambda -\beta x_1{x}_3-\delta x_1-u{x}_1$

${\stackrel{\dot{}}{x}}_2=\beta x_1{x}_3-\epsilon x_2-\delta x_2$

${\stackrel{\dot{}}{x}}_3=\epsilon x_2-\gamma x_3-\delta x_3$

${\stackrel{\dot{}}{x}}_4=\gamma x_3-\delta x_4$

The adjoint variables are:

${\stackrel{\dot{}}{\lambda }}_1={\lambda }_1\left(u+\beta x_3+\delta \right)-{\lambda }_2\beta x_3$

${\stackrel{\dot{}}{\lambda }}_2={\lambda }_2\left(\epsilon +\delta \right)-{\lambda }_3\epsilon$

${\stackrel{\dot{}}{\lambda }}_3=-A-{\lambda }_4\gamma +{\lambda }_3\left(\gamma +\delta \right)-{\lambda }_2\beta x_1+{\lambda }_1\beta x_1$

${\stackrel{\dot{}}{\lambda }}_4={\lambda }_4\delta$

with transversality conditions $\lambda \left(t_f\right)=0$. Figure 8 displays the optimal curves for the states variables and optimal control corresponding to the S-E-I-R model (5).

Figure 8. Solutions of optimal control problem for S-E-I-R model (3.1). u := Vaccination related variable, S := Susceptible Population, E := Exposed Population, I := Infective Population, R := Removed Population.

7. Global Sensitivity Analysis

Global sensitivity analysis allows us to change all parameters simultaneously over the entire parameter interval. This is a way to evaluate the relative effects of each input parameter and also to identify the interactions between parameters to the model output. In global sensitivity analysis we determine that with variation of input parameters in a certain range, which parameters and interactions have the most influential impact on the overall behavior of our model [20] - [26].

There are several types of global sensitivity analyses, such as weighted average of local sensitivity analysis, partial rank correlation coefficient, multi parametric sensitivity analysis, Fourier amplitude sensitivity analysis (FAST) and Sobol’s method, which can be used for systems pharmacology models [20]. The Latin hypercube sampling (LHS) method has been used frequently for global sensitivity analysis. There are also some other methods for calculating main effect and total effect sensitivity indices and one of the most important one is the method of Sobol [25].

LHS method is a sampling method and requires fewer samples compare to simple random sampling to achieve the same accuracy [20]. In LHS method, we divide the random parameter distributions into N equal probability intervals. Here, N is the sample size. The choice for N should be at least k + 1, where k is the number of parameters which are varied. For the case that the interval of variation for some parameter is very large, the sampling can be done on a log form.

In LHS method, sampling is independent for each parameter and can be done by randomly selecting values from each pdf. We may sample each interval once for each parameter without any replacement. The LHS matrix is consisting of N rows corresponding to the number of simulations or sample size and also it includes k columns corresponding to the number of varied parameters. Then, N model solutions may be simulated, using each combination of parameter values which they represent each row of the LHS matrix [20].

7.1. Partial Rank Correlation Coefficient (PRCC) Results for S-I-R Model (2.1)

Here, a parameter sensitivity analysis has being conducted to identify the biological parameters that have the most significant effect on our model system by the LHS Monte Carlo method using PRCC with uniform distributions for the 95 percent confidence intervals. The global sensitivity results with p-values corresponding to S compartment, I compartment and R compartment have been demonstrated in Figure 9.

7.2. Partial Rank Correlation Coefficient (PRCC) Results for S-E-I-R Model (3.1)

According to LHS, we simulated the responses of the model for each organ by randomly selecting values for the parameter set from the 95 percent confidence intervals. These analyses were done by developing a LHS/PRCC method with uniform distributions for the 95 percent confidence intervals. We found that some parameters illustrate significant performance in terms of sensitivity of the output to the variations of these parameters in some compartments while they do not have this effect for others. These results have been depicted in Figure 10, are statistically significant with p-values much smaller than 0.01.

8. Conclusions

Infectious diseases can be defined as diseases that can be transmitted from human to human, from human to animal, or from animal to animal. The mathematical modeling of infectious disease spread has been studied for many years

Figure 9. Global uncertainty and sensitivity analysis of calculated different parameters for S-I-R model (2.1).

Figure 10. Global uncertainty and sensitivity analysis of calculated different parameters for S-E-I-R model (3.1).

and recently it has been widely discussed due to the spread of the COVID-19 pandemic. To build up an appropriate infectious disease dynamic model, we may need to use a system of ordinary differential equations that cover the spread process, spread law, and spread trend of infectious diseases.

In this paper, we considered the S-I-R and S-E-I-R models and for these, we could develop some analytical results which can be useful in studying the simple epidemics. We displayed the evolution of these two compartmental models over time, Susceptible-Infected-Recovered and Susceptible-Exposed-Infected-Recovered for interesting values of parameters. We followed the optimal control perspective to study these models and because of the complexity of the presented optimal control problems, we could no longer solve them analytically and we ended up looking at the numerical solutions. The optimal curves for the states variables and optimal control were obtained and demonstrated for each control problem separately.

An uncertainty analysis can be applied on the epidemiological models to investigate the uncertainty in system output that is generated from uncertainty in parameter inputs. Sensitivity analysis assesses how variations in model outputs can be apportioned, qualitatively or quantitatively, to different inputs. The final objective of this study was to determine the key parameters in spread of infectious diseases using sampling-based method (Partial Rank Correlation Coefficient-PRCC). In this research, we applied LHS/PRCC method with uniform distributions for the 95 percent confidence intervals on the model Equation (2.1) and Equation (3.1). As we have seen, some parameters have positively and some others negatively affected the spread of disease.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Urabe, C.T., Tanaka, G., Aihara, K. and Mimura, M. (2017) Parameter Scaling for Epidemic Size in a Spatial Epidemic Model with Mobile Individuals. PLoS ONE, 11, e0168127. https://doi.org/10.1371/journal.pone.0168127
[2]	Li, G.-H. and Zhang, Y.-X. (2016) Dynamic Behaviors of a Modified SIR Model in Epidemic Diseases Using Nonlinear Incidence and Recovery Rates. PLoS ONE, 12, e0175789. https://doi.org/10.1371/journal.pone.0175789
[3]	Huang, C.X., Cao, J., Wen, F.H. and Yang, X.G. (2016) Stability Analysis of SIR Model with Distributed Delay on Complex Networks. PLoS ONE, 11, e0158813. https://doi.org/10.1371/journal.pone.0158813
[4]	Bauch, C., d’Onofrio, A. and Manfredi, P. (2013) Behavioral Epidemiology of Infectious Diseases: An Overview. In: Modeling the Interplay between Human Behavior and the Spread of Infectious Diseases, Springer, Berlin, 1-19. https://doi.org/10.1007/978-1-4614-5474-8_1
[5]	Kirk, D.E. (2004) Optimal Control Theory: An Introduction. Courier Corporation, North Chelmsford.
[6]	Bryson, A.E. (1996) Optimal Control-1950 to 1985. IEEE Control Systems Magazine, 16, 26-33. https://doi.org/10.1109/37.506395
[7]	Rodrigues, H.S. (2014) Optimal Control and Numerical Optimization Applied to Epidemiological Models.
[8]	Pontryagin, L.S. (2018) Mathematical Theory of Optimal Processes. Routledge, Abingdon-on-Thames. https://doi.org/10.1201/9780203749319
[9]	Leitmann, G. (2013) The Calculus of Variations and Optimal Control: An Introduction. Springer Science & Business Media, Berlin, 24.
[10]	Rao, A.V. (2009) A Survey of Numerical Methods for Optimal Control. Advances in the Astronautical Sciences, 135, 497-528.
[11]	Betts, J.T. (2010) Practical Methods for Optimal Control and Estimation Using Nonlinear Programming. SIAM, University City. https://doi.org/10.1137/1.9780898718577
[12]	Kermack, W.O. and McKendrick, A.G. (1927) A Contribution to the Mathematical Theory of Epidemics. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, the Royal Society London, 115, 700-721. https://doi.org/10.1098/rspa.1927.0118
[13]	Gallina, R. (2012) Dynamic Models for the Analysis of Epidemic Spreads Modelli dinamici per l’analisi di diffusioni epidemiche.
[14]	Lenhart, S. and Workman, J.T. (2007) Optimal Control Applied to Biological Models. CRC Press, Boca Raton. https://doi.org/10.1201/9781420011418
[15]	Chachuat, B. (2007) Nonlinear and Dynamic Optimization: From Theory to Practice.
[16]	Zabczyk, J. (2008) Mathematical Control Theory, Modern Birkhäuser Classics. Birkhäuser Boston Inc., Boston.
[17]	Clarke, F.H. (1990) Optimization and Nonsmooth Analysis. SIAM, University City. https://doi.org/10.1137/1.9781611971309
[18]	Fleming, W.H., Rishel, R.W., Marchuk, G.I., Balakrishnan, A.V., Borovkov, A.A., Makarov, V.L., Rubinov, A.M., Liptser, R.S., Shiryayev, A.N., Krassovsky, N.N., et al. (1975) Applications of Mathematics. Deterministic and Stochastic Optimal Control. Springer-Verlag, Berlin. https://doi.org/10.1007/978-1-4612-6380-7
[19]	Kamien, M.I. and Schwartz, N.L. (2012) Dynamic Optimization: The Calculus of Variations and Optimal Control in Economics and Management. Courier Corporation, North Chelmsford.
[20]	Marino, S., Hogue, I.B., Ray, C.J. and Kirschner, D.E. (2008) A Methodology for Performing Global Uncertainty and Sensitivity Analysis in Systems Biology. Journal of Theoretical Biology, 254, 178-196. https://doi.org/10.1016/j.jtbi.2008.04.011
[21]	Blower, S.M. and Dowlatabadi, H. (1994) Sensitivity and Uncertainty Analysis of Complex Models of Disease Transmission: An HIV Model, as an Example. International Statistical Review, 3, 229-243. https://doi.org/10.2307/1403510
[22]	Zi, Z.K. (2011) Sensitivity Analysis Approaches Applied to Systems Biology Models. IET Systems Biology, 5, 336-346. https://doi.org/10.1049/iet-syb.2011.0015
[23]	Azizi, T. and Mugabi, R. (2020) Global Sensitivity Analysis in Physiological Systems. Applied Mathematics, 11, 119-136. https://doi.org/10.4236/am.2020.113011
[24]	Dalberg, J., Gimenez, H., Keeley, A., Azizi, T., Xi, X. and Jaberi-Douraki, M. (2019) Local and Global Dynamics of Discrete Type 1 Diabetes Model.
[25]	Sobol, I.M. (1993) Sensitivity Estimates for Nonlinear Mathematical Models. Mathematical Modelling and Computational Experiments, 4, 407-414.
[26]	Saltelli, A., Tarantola, S., Campolongo, F. and Ratto, M. (2004) Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models. John Wiley & Sons, Ltd., Hoboken.