∎

11institutetext: Eleni Zavrakli 22institutetext: Hamilton Institute, Maynooth
University, Co. Kildare, Ireland
I-Form Advanced Manufacturing Research Centre, Ireland
Department of Mathematics and Statistics, Maynooth University, Co. Kildare,
Ireland
22email<EMAIL_ADDRESS>33institutetext: Andrew Parnell 44institutetext:
Hamilton Institute, Maynooth University, Co. Kildare, Ireland
I-Form Advanced Manufacturing Research Centre, Ireland
Department of Mathematics and Statistics, Maynooth University, Co. Kildare,
Ireland 55institutetext: David Malone 66institutetext: Hamilton Institute,
Maynooth University, Co. Kildare, Ireland
Department of Mathematics and Statistics, Maynooth University, Co. Kildare,
Ireland 77institutetext: Ken Duffy 88institutetext: Hamilton Institute,
Maynooth University, Co. Kildare, Ireland 99institutetext: Subhrakanti Dey
1010institutetext: Department of Electrical Engineering, Uppsala University,
Sweden

# Optimal age-specific vaccination control for COVID-19: an Irish case study

Eleni Zavrakli  Andrew Parnell  David Malone  Ken Duffy  Subhrakanti Dey

###### Abstract

The outbreak of a novel coronavirus causing severe acute respiratory syndrome
in December 2019 has escalated into a worldwide pandemic. In this work, we
propose a compartmental model to describe the dynamics of transmission of
infection and use it to obtain the optimal vaccination control. The model
accounts for the various stages of the vaccination and the optimisation is
focused on minimising the infections to protect the population and relieve the
healthcare system. As a case study we selected the Republic of Ireland. We use
data provided by Ireland’s COVID-19 Data-Hub and simulate the evolution of the
pandemic with and without the vaccination in place for two different
scenarios, one representative of a national lockdown situation and the other
indicating looser restrictions in place. One of the main findings of our work
is that the optimal approach would involve a vaccination programme where the
older population is vaccinated in larger numbers earlier while simultaneously
part of the younger population also gets vaccinated to lower the risk of
transmission between groups. We compare our simulated results with that of the
vaccination policy taken by the Irish government to explore the advantages of
our optimisation method. Our comparison suggests that a similar reduction in
cases may have been possible even with a reduced set of vaccinations being
available for use.

###### Keywords:

COVID-19 Vaccination Compartmental model Optimal control

## Declarations

### Funding

This publication has emanated from research supported in part by a research
grant from Science Foundation Ireland (SFI) under Grant Number 16/RC/3872 and
is co-funded under the European Regional Development Fund. This work was also
supported by a Science Foundation Ireland COVID-19 Rapid Response Grant number
20/COV/0081. Andrew Parnell’s work was additionally supported by: a Science
Foundation Ireland Career Development Award (17/CDA/4695); an investigator
award (16/IA/4520); a Marine Research Programme funded by the Irish
Government, co-financed by the European Regional Development Fund (Grant-Aid
Agreement No. PBA/CC/18/01); European Union’s Horizon 2020 research and
innovation programme InnoVar under grant agreement No 818144; SFI Centre for
Research Training in Foundations of Data Science 18CRT/6049, and SFI Research
Centre award 12/RC/2289_P2.

### Conflicts of interest

The authors have no conflicts of interest to declare that are relevant to the
content of this article.

### Code availability

The code used for our simulations was written in R and can be found on:
https://github.com/elenizavrakli/Optimal-Age-Specific-Vaccination-Control

##  1 Introduction

In late 2019, an outbreak of pneumonia of unknown cause was reported in the
city of Wuhan in the Hubei province of China cite-key ; Who1 ;
doi:10.1056/NEJMoa2001017 ; https://doi.org/10.1002/jmv.25678 ; cluster ;
newcovid ; doi:10.1056/NEJMoa2001316 . The virus was named Severe Acute
Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) by the World Health
Organisation Who2  and the disease that it causes is referred to as COVID-19
Who3 . The disease quickly became a source of international worry as it spread
around China with most countries around the world following geogr . By March
2020 most countries in the world had confirmed cases of COVID-19, including
the Republic of Ireland, where the first confirmed case was on February 29th
2020 times , the same day that WHO raised the risk warning for the virus to
”very high” times2 . In the absence of an effective treatment or a vaccine,
governments worldwide started implementing protective measures with most of
them announcing a national lockdown to try and control the spread of the virus
and reduce the strain on their healthcare systems. In Ireland, the first
restrictive measures and social distancing guidelines were first announced on
March 12th 2020 and were initially intended to last until March 29th. While
countries tried to control the virus and protect their citizens, the
scientific community committed to coming up with an effective vaccine to put
an end to the pandemic. At the same time, researchers dedicated themselves to
study, model and predict the evolution of the pandemic as well as investigate
non-pharmaceutical intervention methods to control the spread of the virus
Peng2020.02.16.20023465 ; early-transmission ; zhao2020preliminary ;
wu2020nowcasting ; Liu2020.01.25.919787 ; Shen2020.01.23.916726 ;
kucharski2020early ; data-driven ; godio2020seir ; GARBA2020108441 ;
MANEVSKI2020108466 ; NARAIGH2020108496 ; NGONGHALA2020108364 .

On November 9th 2020 Pfizer and BioNTech pfizer ; pfizer-publ  announced a
vaccine candidate that successfully completed the clinical trials and is 90%
effective in preventing infection from the virus. Shortly after, two more
vaccines were announced, namely the one by Moderna moderna  and the one by
AstraZeneca astra . All three vaccines got approved by the European Union ema-
pfizer ; ema-moderna ; ema-astra  with a fourth vaccine getting granted a
conditional marketing authorisation in March 2021, namely the Jcovden one
(previously known as Janssen) ema-johnson . Worldwide, certain states have
made different decisions regarding available vaccines, as there are more
vaccines available. Since the start of 2021, a vaccination roll-out commenced
in most countries, including Ireland, taking into account the number of
available vaccines and the level of risk different groups of people are
considered to be in. As of 2022, first world countries have for the most part
completed the first round of vaccinations and are in the process of
administering booster shots. On the other hand, certain countries are expected
to gain access to sufficient vaccines as late as 2023. padma2021covid . Given
this situation, a need that naturally arises is that of a way to determine the
optimal vaccination strategy, especially given that the resources to reduce
the severity of the pandemic are still limited ACUNAZEGARRA2021108614 ;
RAO2021108621 . This idea was the main motivation for our study.

Modelling, predicting and controlling the behaviour of epidemics has been a
widely studied area booklenhart2007optimal ; bailey1975mathematical . A very
prominent example is the pandemic influenza, a virus that caused an outbreak
of severe pneumonia in 2009, commonly known as ”swine-flu” chowell2009severe .
A number of models were developed to evaluate the implementation of mitigation
strategies bowman2011evaluation ; patel2005finding ; ferguson-nature ; Lee2012
, with a great focus on optimising these strategies.

In this work, we introduce a new model that describes the transmission
dynamics of the virus among different groups of the population in a more
complete way compared to the well known SEIR model Peng2020.02.16.20023465 ;
bacaer2011short ; Mkendrick1925 ; anderson1979population  which is commonly
used in the study of epidemics. The novelty of our approach lies in the
introduction of new compartments to the model accounting for stages of the
vaccination process as well as two different age groups. Using this model as
our baseline, we explore a way to determine the optimal vaccination strategy,
applying optimal control theory methods Fleming1975 ; ogata2010modern ;
kirk2004optimal ; BOLTYANSKI1960464 . Using these tools, we studied the
evolution of the pandemic in the Republic of Ireland, starting from January
2021, using data provided by Ireland’s COVID-19 Data-Hub datahub . We obtained
the optimal vaccination strategy based on estimates of the initial conditions,
which is similar to the course of action taken by the state. However, we find
that it is beneficial to start vaccinating people under the age of 65 in
parallel with people in older age groups but in smaller numbers, as opposed to
exclusively vaccinating the older population first. We compare our obtained
strategy with the course of action taken by the state and comment on the
benefits and disadvantages of both approaches.

In Section 2 we introduce our model and express it as a system of ordinary
differential equations. In Section 3 we use optimal control theory techniques
to obtain the optimal vaccination strategy, given our model. In Section 4 we
produce two sets of simulations to compare the evolution of the pandemic with
and without the vaccination strategy in place, under strict and loose
restrictions. In Section 5 we present a comparative study of the strategy
obtained through our optimisation against the strategy chosen by the
government. Finally, in Section 6 we discuss our findings, the possible
drawbacks of our method and extensions to our approach that are worth
exploring.

##  2 Model

The model most commonly used in the study of epidemics was first introduced by
M’Kendrick in Mkendrick1925 ; anderson1979population  and is known as the SIR
model. This model suggests that every member of a population can be considered
to belong in one of the three compartments: Susceptible (S), Infectious (I) or
Recovered (R). In the study of COVID-19, an extra compartment is usually added
to the SIR model, namely the Exposed (E) compartment, consisting of the people
who have been in contact with the virus yet have not developed any symptoms.
This is known as the SEIR model Peng2020.02.16.20023465 . All compartments can
be indexed by time (t), as they are expressing the number of individuals in
each compartment for each moment in time. In a closed population with no
births and no deaths, the model can be expressed as a system of ordinary
differential equations as follows:

$\displaystyle\frac{dS_{t}}{dt}$ $\displaystyle=$ $\displaystyle-\frac{\beta
I_{t}}{N}S_{t}$ $\displaystyle\frac{dE_{t}}{dt}$ $\displaystyle=$
$\displaystyle\frac{\beta I_{t}}{N}S_{t}-\sigma E_{t}$
$\displaystyle\frac{dI_{t}}{dt}$ $\displaystyle=$ $\displaystyle\sigma
E_{t}-\gamma I_{t}$ $\displaystyle\frac{dR_{t}}{dt}$ $\displaystyle=$
$\displaystyle\gamma I_{t}$

where $N=S_{t}+E_{t}+I_{t}+R_{t}$ and it is constant since we are dealing with
a closed population. $\beta$ is the infectious rate, which expresses the
probability that a susceptible person gets infected by an infectious person,
$\sigma$ is the incubation rate, meaning the rate at which an exposed
individual becomes infectious and $\gamma$ is the recovery rate calculated as
the inverse of the average time of infection.

In this study we use a model that takes the idea of the SEIR model one step
further, with additional compartments related to whether an individual has
been vaccinated and whether the vaccination was effective, hence protecting
the individual. The new set of states consists of the following: Susceptible
not yet vaccinated (S); Received the vaccine; waiting for it to take effect
(V); Received vaccine but it was not effective (N); Susceptible, refusing or
unable to receive the vaccine (U); Exposed, infectious but still asymptomatic
(E); Infectious symptomatic (I); Recovered or deceased (R); Protected from
vaccine (P).

A full table of the notations used throughout the paper can be found in A. In
addition, we assume that a susceptible individual may get infected by exposed
as well as infectious individuals. Furthermore, we consider two different age
groups, those over 65 years old and those younger than 65. The reasons for
this distinction can be found in mueller2020does , the main ones being that
80% of the hospitalised individuals are over the age of 65 and this age group
has 23 times higher risk of death compared to those under 65. This results in
two identical models, one for the younger population and one for the older
which allows for the introduction of different parameters to account for the
way the virus affects people of different ages. Figure 1 is a depiction of our
model for an individual age cohort (for brevity we do not show the full model
configuration). The state configuration is exactly the same in both cases with
the difference being in the transition rates from one state to the other.

Figure 1: Vaccination Compartmental model The model compartments (states) are
indicated by the coloured circles. Blue indicates states at which one might
get infected by an exposed or infectious person and green indicates the states
at which a person is considered to be safe from infection. The parameters
describing the transmission rates are indicated on the arcs between the
states. The models for the over and under 65 populations (indicated by O and Y
subscripts) are identical with the only differences being the age specific
control functions $u_{O}(t),u_{Y}(t)$ which express the percentage of the
susceptible population to be vaccinated at each time point, and the rates of
infection from an exposed or infectious individual
$B_{O}=\frac{\beta_{OO}(E_{O}(t)+I_{O}(t))}{T_{O}}+\frac{\beta_{YO}(E_{Y}(t)+I_{Y}(t))}{T_{Y}}$
and
$B_{Y}=\frac{\beta_{OY}(E_{O}(t)+I_{O}(t))}{T_{O}}+\frac{\beta_{YY}(E_{Y}(t)+I_{Y}(t))}{T_{Y}}$.

We distinguish the states of the two models by adding an O or a Y as an
subscript to the state name to indicate the older or younger than 65 age group
respectively e.g. $S_{O}$, $V_{Y}$ etc. Since the states are time dependent,
we can write them in function notation e.g. $S_{O}(t)$, $V_{Y}(t)$, however,
for brevity, we will omit the $(t)$ in the equations that follow. We also use
the notation $T_{O},T_{Y}$ to indicate the total number of people in each age
group which is assumed to not be time dependent. The two populations influence
each other in the sense that a susceptible person in any of the two groups may
be infected by an exposed or infectious person from either group. For brevity,
we will indicate the older than 65 age group as o65 and the younger than 65
age group as y65. We will be using upper case Roman characters to indicate the
model states and lower case Greek characters to indicate model parameters.

The model is expressed as a system of ordinary differential equations (ODEs),
each of them describing the evolution through time of one of the states. The
ODEs describing the dynamics of our model are given by:

$\displaystyle\frac{dS_{O}}{dt}$ $\displaystyle=$
$\displaystyle-\left(\frac{\beta_{OO}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YO}(E_{Y}+I_{Y})}{T_{Y}}\right)S_{O}-u_{O}S_{O}$
$\displaystyle\frac{dS_{Y}}{dt}$ $\displaystyle=$
$\displaystyle-\left(\frac{\beta_{OY}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YY}(E_{Y}+I_{Y})}{T_{Y}}\right)S_{Y}-u_{Y}S_{Y}$
$\displaystyle\frac{dV_{O}}{dt}$ $\displaystyle=$ $\displaystyle
u_{O}S_{O}-\left(\frac{\beta_{OO}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YO}(E_{Y}+I_{Y})}{T_{Y}}\right)V_{O}$
$\displaystyle-(1-\alpha_{V})\gamma_{V}V_{O}-\gamma_{V}\alpha_{V}V_{O}$
$\displaystyle\frac{dV_{Y}}{dt}$ $\displaystyle=$ $\displaystyle
u_{Y}S_{Y}-\left(\frac{\beta_{OY}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YY}(E_{Y}+I_{Y})}{T_{Y}}\right)V_{Y}$
$\displaystyle-(1-\alpha_{V})\gamma_{V}V_{Y}-\gamma_{V}\alpha_{V}V_{Y}$
$\displaystyle\frac{dN_{O}}{dt}$ $\displaystyle=$
$\displaystyle(1-\alpha_{V})\gamma_{V}V_{O}-\left(\frac{\beta_{OO}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YO}(E_{Y}+I_{Y})}{T_{Y}}\right)N_{O}$
$\displaystyle\frac{dN_{Y}}{dt}$ $\displaystyle=$
$\displaystyle(1-\alpha_{V})\gamma_{V}V_{Y}-\left(\frac{\beta_{OY}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YY}(E_{Y}+I_{Y})}{T_{Y}}\right)N_{Y}$
$\displaystyle\frac{dU_{O}}{dt}$ $\displaystyle=$
$\displaystyle-\left(\frac{\beta_{OO}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YO}(E_{Y}+I_{Y})}{T_{Y}}\right)U_{O}$
$\displaystyle\frac{dU_{Y}}{dt}$ $\displaystyle=$
$\displaystyle-\left(\frac{\beta_{OY}\cdot(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YY}\cdot(E_{Y}+I_{Y})}{T_{Y}}\right)U_{Y}$
$\displaystyle\frac{dE_{O}}{dt}$ $\displaystyle=$
$\displaystyle\left(\frac{\beta_{OO}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YO}(E_{Y}+I_{Y})}{T_{Y}}\right)\left(S_{O}+V_{O}+N_{O}+U_{O}\right)$
$\displaystyle-\gamma_{E}E_{O}$ $\displaystyle\frac{dE_{Y}}{dt}$
$\displaystyle=$
$\displaystyle\left(\frac{\beta_{OY}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YY}(E_{Y}+I_{Y})}{T_{Y}}\right)\left(S_{Y}+V_{Y}+N_{Y}+U_{Y}\right)$
$\displaystyle-\gamma_{E}E_{Y}$ $\displaystyle\frac{dI_{O}}{dt}$
$\displaystyle=$ $\displaystyle\gamma_{E}E_{O}-\gamma_{I}I_{O}$
$\displaystyle\frac{dI_{Y}}{dt}$ $\displaystyle=$
$\displaystyle\gamma_{E}E_{Y}-\gamma_{I}I_{Y}$
$\displaystyle\frac{dR_{O}}{dt}$ $\displaystyle=$
$\displaystyle\gamma_{I}I_{O}$ $\displaystyle\frac{dR_{Y}}{dt}$
$\displaystyle=$ $\displaystyle\gamma_{I}I_{Y}$
$\displaystyle\frac{dP_{O}}{dt}$ $\displaystyle=$
$\displaystyle\gamma_{V}\alpha_{V}\cdot V_{O}$
$\displaystyle\frac{dP_{Y}}{dt}$ $\displaystyle=$
$\displaystyle\gamma_{V}\alpha_{V}\cdot V_{Y}$ (1)

where $T_{O},T_{Y}$ denote the total number of people in the over and under 65
years old populations respectively. These populations are divided in those
willing and those unwilling to get vaccinated ($S_{O},S_{Y},U_{O},U_{Y}$).
This division happens based on estimates of the percentages of people who
refuse the vaccine in each age group. The only option for individuals in the
$U$ compartments, apart from staying there if the pandemic ends early, is to
eventually get infected and be moved to the exposed compartments $E$.

The terms $u_{O},u_{Y}$ are the control functions and represent the percentage
of the old and young population respectively to be vaccinated at each time
point, which means that they are also time dependent. Each control function is
applied to the respective susceptible population ($S_{O},S_{Y}$), moving that
percentage of the population to the respective vaccinated compartment
($V_{O},V_{Y}$). The optimal values for $u_{O},u_{Y}$ are obtained through the
application of the Optimal Control techniques Fleming1975 ;
neustadt1976optimization  discussed in the next section. Any person in the
states $(S,V,U,N)$ is considered susceptible to getting infected by an
infectious or exposed person $(E,I)$. That is because they are either not yet
chosen for vaccination, they received the vaccine and are waiting for it to
take effect, they received the vaccine and it was ineffective, or they chose
not to get vaccinated. The terms $\beta_{ij},i,j\in\\{O,Y\\}$ describe the
transmission rates of the infection between the age groups and they are
analogous to the $\beta$ parameter in the standard SEIR model. Specifically:

  * •

$\beta_{OO}=$ rate at which an o65 susceptible person becomes infected by an
o65 exposed or o65 infected person

  * •

$\beta_{YO}=$ rate at which an o65 susceptible person becomes infected by an
y65 exposed or y65 infected person

  * •

$\beta_{OY}=$ rate at which a y65 susceptible person becomes infected by an
o65 exposed or o65 infected person

  * •

$\beta_{YY}=$ rate at which a y65 susceptible person becomes infected by an
y65 exposed or y65 infected person

There are three rates taken into consideration in the model, the first being
$\gamma_{E}$ which is the rate at which an exposed person becomes symptomatic.
This rate can be calculated as the inverse of the mean holding time to develop
symptoms and become infectious. The next rate is $\gamma_{I}$ which is the
rate at which an infectious person recovers and can be calculated as the
inverse of the mean time to recovery. The final rate considered in the model
is $\gamma_{V}$ which is the rate at which vaccination becomes effective and
can be calculated as the inverse of the mean holding time until protected from
the vaccine. The last parameter influencing the model is $\alpha_{V}$ which is
the vaccine effectiveness. This is expressed as the percentage of people who
become protected from the vaccine, out of those who have received it.

##  3 Optimal control

Optimal control theory Fleming1975 ; neustadt1976optimization ;
kirk2004optimal  is the study of strategies to obtain the control function
that optimises a certain objective over a time horizon, possibly under
suitable constraints. These types of techniques have been widely adopted to
biological systems in general booklenhart2007optimal  but also more
specifically to obtain optimal vaccination strategies when dealing with
viruses and epidemics optimal-epidemics ; Lee2012 ; sharomi2017optimal ;
kirschner1997optimal . More specifically, Pontryagin’s Maximum Principle
pontryagin2018mathematical ; BOLTYANSKI1960464 ; lee1967foundations  is the
main tool that is used when dealing with a problem whose dynamics are
described by a set of Ordinary Differential Equations.

An optimisation problem can be expressed as the problem of minimising an
objective (cost) functional under certain constraints. Let $f(t,x,u)$ denote
the objective functional. Also, let $g(t,x,u)$ denote the state equation (or
set of state equations) of our system. Using $f$ and $g$ we can form a
Hamiltonian function kirk2004optimal ; ogata2010modern  as follows:

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

where $\lambda$ is a continuous function of time$(t)$. For simplicity, we will
write $\lambda(t)$ as just $\lambda$ in the expressions that follow.

For the COVID-19 control problem and our spcific model, we can express the
goal of our optimisation as the minimisation of the number of infectious
individuals, at a minimal cost via vaccination, within a certain time frame
$[0,T]$. That goal can be expressed with the help the following objective
functional to be minimised:

$\mathcal{F}\left(U(t)\right)=\int_{0}^{T}\left[I_{O}(t)+I_{Y}(t)+\frac{W_{O}}{2}u_{O}^{2}(t)+\frac{W_{Y}}{2}u_{Y}^{2}(t)\right]dt$
(2)

where $U(t)=\left(u_{O}(t),u_{Y}(t)\right)$ and $W_{O},W_{Y}\geq 0$ are the
age specific weight constants enforcing the severity of the optimisation
constraints. This type of choice of the objective functional as a non-negative
increasing function of the control inputs, is common practice in optimal
control applications Lee2012 ; sharomi2017optimal ; kirk2004optimal . The
control functions are squared in order to ensure the convexity of the
functional, guaranteeing a sufficient condition for the existence of an
optimal solution (to be proved later in detail). Consider $X(t)=$
$\left(\right.S_{O}(t),$ $V_{O}(t),$ $N_{O}(t),$ $U_{O}(t),$ $E_{O}(t),$
$I_{O}(t),$ $R_{O}(t),$ $P_{O}(t),$ $S_{Y}(t),$ $V_{Y}(t),$ $N_{Y}(t),$
$U_{Y}(t),$ $E_{Y}(t),$ $I_{Y}(t),$ $R_{Y}(t),$ $P_{Y}(t)\left.\right)$. We
are looking for the optimal pair of solutions
$\left(U^{*}(t),X^{*}(t)\right)$, i.e. the optimal control $U^{*}$ and the
corresponding trajectory $X^{*}$ when $U^{*}$ is applied, such that

$\mathcal{F}(U^{*}(t))=\min_{\Omega}\mathcal{F}(U(t))$ (3)

where $\Omega=\left\\{U(t)\in L^{2}(O,T)^{2}\|a\leq u_{O}(t),u_{Y}(t)\leq
b,t\in[0,T]\right\\}$, $a$ and $b$ are the upper and lower bounds for the
control function and can usually be expressed as real values, and $T$ is the
time horizon for our optimisation.

There are a few different approaches that can be taken in defining the
constrained optimisation problem and are worth mentioning. We can include
weighting factors to the infectious compartments ($I_{O},I_{Y}$) in the
objective functional $\mathcal{F}$ that correspond to the mortality rates for
each age group or more factors that generally model the cost that large number
of infections can have in the healthcare system. Additionally, when it comes
to bounding the control functions, we can chose a constraint of the form
$a\leq u_{O}(t)+u_{Y}(t)\leq b$ which can be interpreted as the total
percentage of vaccinations being bounded as opposed to bounding them per age
group. This kind of constraint would result in an extra term in the objective
functional and the Hamiltonian function having the form
$\lambda\cdot(u_{O}+u_{Y})$.

The following theorem describes the main conditions that when satisfied can
lead us to the optimal control that solves the minimisation problem.

######  Theorem 3.1 (Pontryagin’s maximum principle)

Given the Hamiltonian

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

then the following conditions are satisfied by the optimal control
$u^{*}\in\Omega=\left\\{U(t)\in L^{2}(O,T)^{2}\|a\leq u_{O}(t),u_{Y}(t)\leq
b,t\in[0,T]\right\\}$:

$\displaystyle\frac{\partial H}{\partial u}\left\\{\begin{matrix}=&0&\text{ at
}u^{*}\text{ if }&a<u^{*}<b\\\ \geq&0&\text{ at }u^{*}\text{ if }&u^{*}=a\\\
\leq&0&\text{ at }u^{*}\text{ if }&u^{*}=b\end{matrix}\right.$ Optimality
Condition $\displaystyle\dot{\lambda}=-\frac{\partial H}{\partial x}$ Adjoint
Equation $\displaystyle\lambda(T)=0$ Transversality Condition
$\displaystyle\dot{x}=g(t,x,u),$ $\displaystyle x(0)=x_{0}$ Dynamics of state
equation

where $\dot{\lambda}=\partial\lambda/\partial t$ and $\dot{x}=\partial
x/\partial t$.

In order to apply the maximum principle we first need to define the
Hamiltonian function, omitting for brevity the time dependence from the state
notations:

$\displaystyle H$ $\displaystyle=$
$\displaystyle\left[I_{O}+I_{Y}+\frac{W_{O}}{2}u_{O}^{2}+\frac{W_{Y}}{2}u_{Y}^{2}\right]$
$\displaystyle+$
$\displaystyle\lambda_{S_{O}}\left\\{-\left(\frac{\beta_{OO}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YO}(E_{Y}+I_{Y})}{T_{Y}}\right)S_{O}-u_{O}S_{O}\right\\}$
$\displaystyle+$
$\displaystyle\lambda_{S_{Y}}\left\\{-\left(\frac{\beta_{OY}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YY}(E_{Y}+I_{Y})}{T_{Y}}\right)S_{Y}-u_{Y}S_{Y}\right\\}$
$\displaystyle+$
$\displaystyle\lambda_{V_{O}}\left\\{u_{O}S_{O}-\left(\frac{\beta_{OO}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YO}(E_{Y}+I_{Y})}{T_{Y}}\right)V_{O}-\gamma_{V}V_{O}\right.$
$\displaystyle\left.-\alpha_{V}V_{O}\right\\}$ $\displaystyle+$
$\displaystyle\lambda_{V_{Y}}\left\\{u_{Y}S_{Y}-\left(\frac{\beta_{OY}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YY}(E_{Y}+I_{Y})}{T_{Y}}\right)V_{Y}-\gamma_{V}V_{Y}\right.$
$\displaystyle\left.-\alpha_{V}V_{Y}\right\\}$ $\displaystyle+$
$\displaystyle\lambda_{N_{O}}\left\\{\gamma_{V}V_{O}-\left(\frac{\beta_{OO}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YO}(E_{Y}+I_{Y})}{T_{Y}}\right)N_{O}\right\\}$
$\displaystyle+$
$\displaystyle\lambda_{N_{Y}}\left\\{\gamma_{V}V_{Y}-\left(\frac{\beta_{OY}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YY}(E_{Y}+I_{Y})}{T_{Y}}\right)N_{Y}\right\\}$
$\displaystyle-$
$\displaystyle\lambda_{U_{O}}\left(\frac{\beta_{OO}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YO}(E_{Y}+I_{Y})}{T_{Y}}\right)U_{O}$
$\displaystyle-$
$\displaystyle\lambda_{U_{Y}}\left(\frac{\beta_{OY}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YY}(E_{Y}+I_{Y})}{T_{Y}}\right)U_{Y}$
$\displaystyle+$
$\displaystyle\lambda_{E_{O}}\left\\{\left(\frac{\beta_{OO}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YO}(E_{Y}+I_{Y})}{T_{Y}}\right)\left(S_{O}+V_{O}+N_{O}+U_{O}\right)\right.$
$\displaystyle\left.-\gamma_{E}E_{O}\right\\}$ $\displaystyle+$
$\displaystyle\lambda_{E_{Y}}\left\\{\left(\frac{\beta_{OY}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YY}(E_{Y}+I_{Y})}{T_{Y}}\right)\left(S_{Y}+V_{Y}+N_{Y}+U_{Y}\right)\right.$
$\displaystyle\left.-\gamma_{E}E_{Y}\right\\}$ $\displaystyle+$
$\displaystyle\lambda_{I_{O}}\left\\{\gamma_{E}E_{O}-\gamma_{I}I_{O}\right\\}$
$\displaystyle+$
$\displaystyle\lambda_{I_{Y}}\left\\{\gamma_{E}E_{Y}-\gamma_{I}I_{Y}\right\\}$

where $\Lambda:=$ $\left\\{\right.\lambda_{S_{O}},$ $\lambda_{S_{Y}},$
$\lambda_{V_{O}},$ $\lambda_{V_{Y}},$ $\lambda_{N_{O}},$ $\lambda_{N_{Y}},$
$\lambda_{U_{O}},$ $\lambda_{U_{Y}},$ $\lambda_{E_{O}},$ $\lambda_{E_{Y}},$
$\lambda_{I_{O}},$ $\lambda_{I_{Y}}\left.\right\\}$ is a set of continuous
functions of time ($\Lambda(t)$), given the states and controls. These
functions have a key role in our optimisation technique.

Applying each of the conditions of the maximum principle (3.1) we obtain the
following:

  * •

From the optimality condition, the optimal control strategies are given by:

$u_{O}^{*}(t)=\min\left\\{\max\left\\{a,\frac{S_{O}(t)}{W_{O}}\left(\lambda_{S_{O}}(t)-\lambda_{V_{O}}(t)\right)\right\\},b\right\\}$
(4)

$u_{Y}^{*}(t)=\min\left\\{\max\left\\{a,\frac{S_{Y}(t)}{W_{Y}}\left(\lambda_{S_{Y}}(t)-\lambda_{V_{Y}}(t)\right)\right\\},b\right\\}$
(5)

  * •

Adjoint Equation: $\dot{\Lambda}=-\partial H/\partial X$ which results in the
following system of ODEs:

$\displaystyle\dot{\lambda}_{S_{O}}$ $\displaystyle=$
$\displaystyle\left(\frac{\beta_{OO}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YO}(E_{Y}+I_{Y})}{T_{Y}}+u_{O}\right)\lambda_{S_{O}}-u_{O}\lambda_{V_{O}}$
$\displaystyle-$
$\displaystyle\left(\frac{\beta_{OO}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YO}(E_{Y}+I_{Y})}{T_{Y}}\right)\lambda_{E_{O}}$
$\displaystyle\dot{\lambda}_{S_{Y}}$ $\displaystyle=$
$\displaystyle\left(\frac{\beta_{OY}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YY}(E_{Y}+I_{Y})}{T_{Y}}+u_{Y}\right)\lambda_{S_{Y}}-u_{Y}\lambda_{V_{Y}}$
$\displaystyle-$
$\displaystyle\left(\frac{\beta_{OY}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YY}(E_{Y}+I_{Y})}{T_{Y}}\right)\lambda_{E_{Y}}$
$\displaystyle\dot{\lambda}_{V_{O}}$ $\displaystyle=$
$\displaystyle\left(\frac{\beta_{OO}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YO}(E_{Y}+I_{Y})}{T_{Y}}+\gamma_{V}\right)\lambda_{V_{O}}-(1-\alpha_{V})\gamma_{V}\lambda_{N_{O}}$
$\displaystyle-$
$\displaystyle\left(\frac{\beta_{OO}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YO}(E_{Y}+I_{Y})}{T_{Y}}\right)\lambda_{E_{O}}$
$\displaystyle\dot{\lambda}_{V_{Y}}$ $\displaystyle=$
$\displaystyle\left(\frac{\beta_{OY}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YY}(E_{Y}+I_{Y})}{T_{Y}}+\gamma_{V}\right)\lambda_{V_{Y}}-(1-\alpha_{V})\gamma_{V}\lambda_{N_{Y}}$
$\displaystyle-$
$\displaystyle\left(\frac{\beta_{OY}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YY}(E_{Y}+I_{Y})}{T_{Y}}\right)\lambda_{E_{Y}}$
$\displaystyle\dot{\lambda}_{N_{O}}$ $\displaystyle=$
$\displaystyle\left(\frac{\beta_{OO}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YO}(E_{Y}+I_{Y})}{T_{Y}}\right)(\lambda_{N_{O}}-\lambda_{E_{O}})$
$\displaystyle\dot{\lambda}_{N_{Y}}$ $\displaystyle=$
$\displaystyle\left(\frac{\beta_{OY}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YY}(E_{Y}+I_{Y})}{T_{Y}}\right)(\lambda_{N_{Y}}-\lambda_{E_{Y}})$
$\displaystyle\dot{\lambda}_{U_{O}}$ $\displaystyle=$
$\displaystyle\left(\frac{\beta_{OO}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YO}(E_{Y}+I_{Y})}{T_{Y}}\right)(\lambda_{U_{O}}-\lambda_{E_{O}})$
$\displaystyle\dot{\lambda}_{U_{Y}}$ $\displaystyle=$
$\displaystyle\left(\frac{\beta_{OY}(E_{O}+I_{O})}{T_{O}}+\frac{\beta_{YY}(E_{Y}+I_{Y})}{T_{Y}}\right)(\lambda_{U_{Y}}-\lambda_{E_{Y}})$
$\displaystyle\dot{\lambda}_{E_{O}}$ $\displaystyle=$
$\displaystyle\frac{\beta_{OO}S_{O}}{T_{O}}\left(\lambda_{S_{O}}-\lambda_{E_{O}}\right)+\frac{\beta_{OO}V_{O}}{T_{O}}\left(\lambda_{V_{O}}-\lambda_{E_{O}}\right)$
$\displaystyle+$
$\displaystyle\frac{\beta_{OO}N_{O}}{T_{O}}\left(\lambda_{N_{O}}-\lambda_{E_{O}}\right)+\frac{\beta_{OO}U_{O}}{T_{O}}\left(\lambda_{U_{O}}-\lambda_{E_{O}}\right)$
$\displaystyle+$
$\displaystyle\frac{\beta_{OY}S_{Y}}{T_{O}}\left(\lambda_{S_{Y}}-\lambda_{E_{Y}}\right)+\frac{\beta_{OY}YV_{t}}{T_{O}}\left(\lambda_{V_{Y}}-\lambda_{E_{Y}}\right)$
$\displaystyle+$
$\displaystyle\frac{\beta_{OY}N_{Y}}{T_{O}}\left(\lambda_{N_{Y}}-\lambda_{E_{Y}}\right)+\frac{\beta_{OY}U_{Y}}{T_{O}}\left(\lambda_{U_{Y}}-\lambda_{E_{Y}}\right)$
$\displaystyle+$
$\displaystyle\gamma_{E}\lambda_{E_{O}}-\gamma_{E}\lambda_{I_{O}}$
$\displaystyle\dot{\lambda}_{E_{Y}}$ $\displaystyle=$
$\displaystyle\frac{\beta_{YO}S_{O}}{T_{Y}}\left(\lambda_{S_{O}}-\lambda_{E_{O}}\right)+\frac{\beta_{YO}V_{O}}{T_{Y}}\left(\lambda_{V_{O}}-\lambda_{E_{O}}\right)$
$\displaystyle+$
$\displaystyle\frac{\beta_{YO}N_{O}}{T_{Y}}\left(\lambda_{N_{O}}-\lambda_{E_{O}}\right)+\frac{\beta_{YO}U_{O}}{T_{Y}}\left(\lambda_{U_{O}}-\lambda_{E_{O}}\right)$
$\displaystyle+$
$\displaystyle\frac{\beta_{YY}S_{Y}}{T_{Y}}\left(\lambda_{S_{Y}}-\lambda_{E_{Y}}\right)+\frac{\beta_{YY}V_{Y}}{T_{Y}}\left(\lambda_{V_{Y}}-\lambda_{E_{Y}}\right)$
$\displaystyle+$
$\displaystyle\frac{\beta_{YY}N_{Y}}{T_{Y}}\left(\lambda_{N_{Y}}-\lambda_{E_{Y}}\right)+\frac{\beta_{YY}U_{Y}}{T_{Y}}\left(\lambda_{U_{Y}}-\lambda_{E_{Y}}\right)$
$\displaystyle+$
$\displaystyle\gamma_{E}\lambda_{E_{Y}}-\gamma_{E}\lambda_{I_{Y}}$
$\displaystyle\dot{\lambda}_{I_{O}}$ $\displaystyle=$
$\displaystyle\frac{\beta_{OO}S_{O}}{T_{O}}\left(\lambda_{S_{O}}-\lambda_{E_{O}}\right)+\frac{\beta_{OO}V_{O}}{T_{O}}\left(\lambda_{V_{O}}-\lambda_{E_{O}}\right)$
$\displaystyle+$
$\displaystyle\frac{\beta_{OO}N_{O}}{T_{O}}\left(\lambda_{N_{O}}-\lambda_{E_{O}}\right)+\frac{\beta_{OO}U_{O}}{T_{O}}\left(\lambda_{U_{O}}-\lambda_{E_{O}}\right)$
$\displaystyle+$
$\displaystyle\frac{\beta_{OY}S_{Y}}{T_{Y}}\left(\lambda_{S_{Y}}-\lambda_{E_{Y}}\right)+\frac{\beta_{OY}V_{Y}}{T_{O}}\left(\lambda_{V_{Y}}-\lambda_{E_{Y}}\right)$
$\displaystyle+$
$\displaystyle\frac{\beta_{OY}N_{Y}}{T_{O}}\left(\lambda_{N_{Y}}-\lambda_{E_{Y}}\right)+\frac{\beta_{OY}U_{Y}}{T_{O}}\left(\lambda_{U_{Y}}-\lambda_{E_{Y}}\right)$
$\displaystyle+$ $\displaystyle\gamma_{I}\lambda_{I_{O}}-1$
$\displaystyle\dot{\lambda}_{I_{Y}}$ $\displaystyle=$
$\displaystyle\frac{\beta_{YO}S_{O}}{T_{Y}}\left(\lambda_{S_{O}}-\lambda_{E_{O}}\right)+\frac{\beta_{YO}V_{O}}{T_{Y}}\left(\lambda_{V_{O}}-\lambda_{E_{O}}\right)$
$\displaystyle+$
$\displaystyle\frac{\beta_{YO}N_{O}}{T_{Y}}\left(\lambda_{N_{O}}-\lambda_{E_{O}}\right)+\frac{\beta_{YO}U_{O}}{T_{Y}}\left(\lambda_{U_{O}}-\lambda_{E_{O}}\right)$
$\displaystyle+$
$\displaystyle\frac{\beta_{YY}S_{Y}}{T_{Y}}\left(\lambda_{S_{Y}}-\lambda_{E_{Y}}\right)+\frac{\beta_{YY}V_{Y}}{T_{Y}}\left(\lambda_{V_{Y}}-\lambda_{E_{Y}}\right)$
$\displaystyle+$
$\displaystyle\frac{\beta_{YY}N_{Y}}{T_{Y}}\left(\lambda_{N_{Y}}-\lambda_{E_{Y}}\right)+\frac{\beta_{YY}U_{Y}}{T_{Y}}\left(\lambda_{U_{Y}}-\lambda_{E_{Y}}\right)$
$\displaystyle+$ $\displaystyle\gamma_{I}\lambda_{I_{Y}}-1$

where $\dot{\lambda}$ denotes $\partial\lambda/\partial t$

  * •

The transversality condition: $\Lambda(T)=0$ gives us a terminal condition for
the system of adjoint ODEs. This condition implies that the adjoint system
should be solved backwards butcher2008numerical  as opposed to the state
system.

  * •

The dynamics of the system are described by the system of state equations (2).

To ensure the existence of a solution to our optimisation problem, we need to
examine the sufficient conditions as stated in Fleming1975 , as also used in
gaff2009optimal . The relevant conditions and the corresponding existence
proof for our model can be found in Appendix B.

Having obtained all the necessary optimality conditions and resulting
equations, Algorithm 1 describes the steps that lead to the computation of an
optimal control trajectory.

Result: Optimal control and trajectory $U^{*},X^{*}$

Input: Initial state $x_{0}$, Parameter values, Time horizon, Initial control,
Convergence threshold $\delta$

1 Solve the system of state ODEs forwards in time to determine the state;

2 Solve the system of adjoint ODEs backwards in time;

3 Update the controls using (4) and (5) ;

4 Compute the error term $\frac{\|x_{k}-x_{k-1}\|_{1}}{\|x_{k-1}\|_{1}}$ ;

If the error term is smaller than some predetermined threshold $\delta$ then
extract the optimal control and optimal trajectory. Otherwise repeat the
process until convergence.

Algorithm 1 Determining the optimal control

##  4 Simulation Setup and Results

We used our model to estimate the evolution of the virus in the Republic of
Ireland, starting from January 2021 using data provided by datahub  regarding
the population of the country, the number of exposed, infected and recovered
people in both age groups. Additionally, we applied the optimal control
strategy on the same data to obtain the effect that the vaccination would have
on the evolution of the pandemic, were it adopted by the state.

Table 1 includes the initial values of the states
$E_{O},I_{O},R_{O},E_{Y},I_{Y},R_{Y}$. The rest of the states can be
calculated with the use of this information and some of the parameter values.
The total susceptible people in each age group are the people not yet
vaccinated along with the people refusing or unable to receive the vaccine
($\textbf{S}:=S+U$). This total of susceptibles for each age group can be
calculated by: $\textbf{S}=T-(E+I+R)$. The separation between not yet
vaccinated ($S$) and unwilling to get vaccinated ($U$) can be expressed as a
percentage for each population, representing the refusal rate for the vaccine.
Let $r_{O}$ and $r_{Y}$ be these percentages in the over 65 and under 65
populations respectively. Then $S_{Y}=(1-r_{Y})\textbf{S}$,
$U_{Y}=r_{Y}\textbf{S}$, $S_{O}=(1-r_{O})\textbf{S}$, $U_{O}=r_{O}\textbf{S}$.

State (Compartment) | Number of people
---|---
Total y65 ($T_{Y}$) | 4000000
Total o65 ($T_{O}$) | 900000
Y65 exposed ($E_{Y}$) | 2000
O65 exposed ($E_{O}$) | 200
Y65 infected ($I_{Y}$) | 2000
O65 infected ($I_{O}$) | 200
Y65 recovered ($R_{Y}$) | 200000
O65 recovered ($R_{O}$) | 100000
Table 1:  Initial values for the model states

We produce two different simulations, to study the evolution of the pandemic,
following the steps described in Algorithm 1. In the first case we assume
strict measures are in place, meaning that the transmission rates of the virus
among individuals is low. In the second case we assume minimal restrictions
are in place, hence the transmission rates are much higher. The parameters
used in both simulations are listed in Table 2.

Parameter Description | Symbol | Case 1 | Case 2
---|---|---|---
Mean no of o65 infected by an o65 | $R0_{OO}$ | 1.2 | 8
Mean no of y65 infected by a y65 | $R0_{YY}$ | 1.2 | 8
Mean no of y65 infected by an o65 | $R0_{OY}$ | 0.9 | 4
Mean no of o65 infected by a y65 | $R0_{YO}$ | 0.9 | 3
Mean holding time exposed in days | $t_{E}$ | 6.6
Mean holding time infected in days | $t_{I}$ | 7.4
Mean holding time vaccinated in days | $t_{V}$ | 14
Vaccine effectiveness | $\alpha_{V}$ | 0.9
Upper bound for control function | $b$ | 0.3
Percentage of o65 people refusing the vaccine | $r_{O}$ | 0.07
Percentage of y65 people refusing the vaccine | $r_{Y}$ | 0.21
Weight constant relating to o65 group | $W_{O}$ | $10^{11}$
Weight constant relating to y65 group | $W_{Y}$ | $10^{11}$
Table 2: Parameter values for both simulations. The mean numbers of infections
from each infected person within each age group and between age groups vary
for the two simulations, while for the rest of the parameters we used the same
values. Case 1 is representative of a situation where strict measures are in
place and the results of the simulation can be found in subsection 4.1. Case 2
is an example of a situation with minimal restrictions and is explored in
subsection 4.2.

The three parameters $t_{E},t_{I},t_{V}$, describe the mean holding times in
each of the states E, I and V. They are the inverses of the rates we have in
our model (2), namely $\gamma_{E},\gamma_{I},\gamma_{V}$. Also, particularly
interesting are the four $RO$ parameters, that are related to how quickly the
infection is transmitted between the members of the entire population. More
specifically, these parameters are directly related to the terms $\beta_{ij}$
in our model which describe the transmission rates, with the help of the three
mean holding times:

$\displaystyle\beta_{OO}=\frac{R0_{OO}}{t_{E}+t_{I}+t_{V}},\beta_{YO}=\frac{R0_{YO}}{t_{E}+t_{I}+t_{V}},$
$\displaystyle\beta_{OY}=\frac{R0_{OY}}{t_{E}+t_{I}+t_{V}},\beta_{YY}=\frac{R0_{YY}}{t_{E}+t_{I}+t_{V}}$
(6)

Two different sets of $R0$ numbers are used to produce our two simulations and
give an insight to the evolution of the pandemic under different levels of
restrictions.

The values we use for the parameters are estimates specific to COVID-19 (mean
holding times) cereda2020early , or based on estimates used in similar studies
of other infectious diseases such as (weight constants) Lee2012 . The upper
bound for the control function ($b$) is chosen arbitrarily, however it is not
a value that is ever obtained for either control function, thanks to the high
weight constants used in the objective functional that model the cost of the
vaccination and the severity of the constraints. An alternative approach would
be for $b$ to bound the sum of the control functions as opposed to each of
them individually. The weight constants $W_{O},W_{Y}$ are chosen to be equal
because we assume no difference in the cost of providing the vaccination to
people in different age groups. However, the two groups can be weighed
differently, explaining the cost of bringing the vaccination to remote
locations, nursing homes etc. The refusal percentages for both age groups
($r_{O},r_{Y}$) we used are based on a survey survey , the results of which
suggest that $77\%$ of the overall population of Ireland are willing to
receive the vaccine and the the willingness is stronger in the over 65
population, with $93\%$ of the group intending to get vaccinated. We included
the percentage of people unsure about the vaccination ($15\%$ overall) in the
refusal percentage ($6\%$) along with the people who claimed to already be
against receiving the vaccination. Finally, as we mentioned in section 3, it
is possible for one more set of parameters to be included in our model and
more specifically the optimisation function, weighing the number of infectious
individuals by the mortality rate for each age group.

###  4.1 Evolution of pandemic under tight restrictions

We simulate the evolution of the pandemic using the initial state values given
in Table 1 and the parameter values given in Table 2. Specifically, the $R0$
numbers we use are given in the column named Case 1 and are an example of the
way the virus is evolving under a strict restriction policy. In that case, the
average number of infections within each age group is 1.2 and between age
groups 0.9.

(a) Evolution of states without the vaccine

(b) Evolution of states with the vaccine

Figure 2:  Evolution of the pandemic with and without vaccination under strict
measures. In the case where no vaccination policy is in place, there is an
early peak for the infectious individuals in both age groups (red and
turquoise curves). In the second graph the peaks are almost indistinguishable.
This is due to the vaccination reducing both the total and the maximum number
of infections that take place, hence flattening the curve. The number of
susceptible people (pink and blue curves) in both age groups declines fast in
the first graph, due to most of them getting infected, while in the second
case the decline is much slower, as many individuals are getting vaccinated
and thus protected from infection.

(a) Evolution of states without the vaccine

(b) Evolution of states with the vaccine

Figure 3: Evolution of the pandemic with and without vaccination under strict
measures. The graphs are produced in vertical log scale to give a clearer view
of the peaks. In the first graph there are clear peaks in the number of
infectious individuals in both age groups while in the second, the curves are
clearly flattened thanks to the vaccination.

Figure 2 shows the difference between the evolution of the pandemic when there
is no control in place and when the optimal control is applied. Due to the
transmission rates being quite low, the increase in the number of infectious
individuals will not be particularly evident, especially when the vaccination
is applied. For that reason we represented the same information also on the
log scale in Figure 3 to give a clearer picture of the difference an optimal
vaccination strategy makes to the evolution of the pandemic. In the case where
no control is applied there is a peak in the number of infectious individuals
in both age groups and as a result, the number of susceptible people drops
fast with the pandemic ending in a very short period of time, namely less than
250 days, with the greatest part of the population having been infected.
Specifically, 80% of the population over 65 (720,249) and 79% of the
population under 65 (3,159,510) gets infected by the virus until the end of
the pandemic. That means that around 79% of the total population of the
country (3,879,759) will get infected by the virus when no vaccination
strategy is in place.

On the other hand, when the optimal vaccination strategy is applied, the
infectious curve is flattened to the point where it becomes nearly
indistinguishable from the horizontal axis. This means that the number of
people to become infected by the virus is substantially reduced as a result of
the protection provided by the vaccine. In particular, only 11.3% (101,764)
and 5.38% (215,132) of the over and under 65 years old groups respectively
will get infected by the virus, resulting in a total 6.47% of the population
(316,896). This is achieved thanks to the 70.58% (635,187) and 17.08%
(683,044) of the over 65 and under 65 populations respectively receiving the
vaccine and being successfully protected by it. In total, that is 26.9% of the
population (1,318,231) being successfully vaccinated.

The difference that the control function makes to the number of infectious
individuals (in the log scale) can be more clearly seen in Figure 4. For both
age groups, the curve is substantially flattened as a result of the
vaccination. This is a visual representation of the reduction in the total
number of infections and also showcases the fact that less people will be
infected at the same time, causing less stress on the healthcare system.

(a) Evolution of infectious in over 65

(b) Evolution of infectious in under 65

Figure 4: Comparison of baseline curves (without vaccine) with optimal
vaccination strategy curves for the infectious populations in both age groups
under strict measures. The vertical scale is in log units. The vaccination
successfully flattens the curve in both cases, resulting in fewer infections
overall as well as fewer simultaneous infections.

(a) Optimal control functions

(b) Optimal control for over65

(c) Optimal control for under65

Figure 5:  Optimal control functions for both age groups. 5(a) presents both
control functions on a vertical log scale. The o65 group gets vaccinated
faster and in higher percentages, while the vaccination of the y65 group
commences but proceeds slower and continues after the pandemic is under
control.

The optimal vaccination strategy for each age group, i.e., the control
functions that resulted from the optimisation technique described in section 3
can be seen in Figure 5. As expected, close to the start of our simulations, a
high percentage of the susceptible population gets vaccinated to quickly get
the pandemic under control, and then that percentage continuously declines. It
is interesting to note how the o65 population gets vaccinated earlier and with
higher percentages while the vaccination of the younger population is
happening slower, with smaller percentages every day, and goes on after the
pandemic is over as can be seen in Figure 5. This is also evident in the
percentages of the people in each age group that get vaccinated by the end of
the simulation, with 70.58% and 17.08% for the over and under 65 groups
respectively, as described earlier.

It is worth noting that the parameter setup used for this simulation is
indicative of a national lockdown situation, which is an unrealistic scenario
for this length of time. For this reason, we looked into a second parameter
setup, discussed in the next section.

###  4.2 Evolution of pandemic under loose restrictions

We now simulate the evolution of the pandemic in Ireland, using the same
initial conditions for the states as defined in Table 1. We are interested in
investigating the effect a different choice of parameters, specifically
different $R0$ numbers, would have on our results. We are looking into values
that are indicative of an extreme scenario where the measures in place are not
strict or even non existent and each infected individual can infect multiple
individuals before recovering. The column titled Case 2 in Table 2 consists of
values chosen for the parameters, while the rest of the parameters remain the
same as in subsection 4.1.

When no vaccination policy is in place, the curves of infectious individuals
in both age groups present a very high peak due to the high transmission rates
as can be seen in Figures 6 and 7. This peak is indicative of a period of time
where a very large number of people are infected at the same time, a situation
that would cause a huge strain on the country’s healthcare system. The
percentages of people that get infected in the over and under 65 group
respectively are 99.995% (899,954) and 99.993% (3,999,713), resulting in a
total 99.993% (4,899,667) of the population. This means that nearly every
single resident of the country will get infected by the virus in a very short
period of time, while at the same time paralysing the healthcare system.

That curve gets flattened, with a lower peak and greater spread, thanks to the
optimal control functions, i.e. the optimal vaccination strategy. This would
offer some relief to the healthcare system, as not all infections will be
active at the same time. The percentages of people that get infected in the
over and under 65 group respectively are 85.29% (767,614) and 99.1%
(3,965,963), resulting in a total 96.6% of the population (4,733,577). These
percentages, though reduced, are still very high and that is due to the fact
that the chosen transmission rates are so high. For the same reason the
percentages of the people that get successfully protected from the vaccine are
relatively low, namely 14.2% (127,871) and 0.1% (3,861) for the over 65 and
under 65 groups respectively, meaning a total 2.7% of the population (131,732)
will be successfully vaccinated. There is an additional reason to why these
percentages are so low. From the moment an individual receives the vaccine,
there is a period equal to the mean holding time vaccinated ($t_{V}$) during
which the person may still get infected. Due to the high transmission rates, a
lot of vaccinated individuals will get infected, hence not making it to the
protected state.

A side by side presentation of the evolution of the pandemic with and without
vaccination is given in Figure 6 and a comparison of the optimised with the
baseline curves for the infectious populations in each age group are given in
Figure 7.

(a) Evolution of states without the vaccine

(b) Evolution of states with the vaccine

Figure 6: Evolution of the pandemic with an without vaccination under loose
measures. Without a vaccination policy in place, the curve of infectious
individuals presents an earlier and higher peak as opposed to the case where
the vaccination strategy is applied and the curve is flattened.

(a) Evolution of infectious in over 65

(b) Evolution of infectious in under 65

Figure 7: Comparison of baseline curves (without vaccine) with optimal
vaccination strategy curves for the infectious populations in both age groups
under loose restrictions.

Figure 8 shows the optimal control functions that can be obtained as the
solution to our optimisation problem. Similar to the previous simulations, the
control function curves present a peak at first, indicating a high percentage
of the population being chosen for vaccination early on and that percentage
getting reduced as time progresses.

(a) Optimal control functions

(b) Optimal control for over65

(c) Optimal control for under65

Figure 8: Optimal control functions for both age groups. 8(a) presents both
control functions in log scale.

##  5 Comparative Study

We now look at the effect that the optimal control obtained from our model
would have on the evolution of the pandemic in Ireland by comparing our
approach to the real infection numbers, given the vaccination strategy taken
by the Irish government. In order for the comparison to be sensible, we limit
our simulations to a time frame of 100 days, thus ensuring that the
transmission rates remain somewhat constant. This is needed because our model
does not account for varying $R0$ numbers, while in reality these numbers can
fluctuate depending on the restrictive measures in place (e.g. lockdowns) and
the preventative measures taken, namely vaccinations. We obtained our data
from datahub , specifically the daily infection and vaccination numbers
starting from January 1st 2021\. These values are shown in Figure9.

(a) Infectious individuals

(b) Vaccines administered

Figure 9: Infection and vaccination numbers for the Republic of Ireland
starting from January 1st 2021

It is evident that there is a substantial peak in the infection numbers during
the first few days which can be attributed to the looser measures in place
before the Christmas holidays of 2020. Considering this, we produce a second
set of comparative plots, starting the simulation just after the peak in
infections, namely on January 18th instead of the 1st. The daily infection and
vaccination numbers for a period of 100 days starting from January 18th can be
seen in Figure10.

(a) Infectious individuals

(b) Vaccines administered

Figure 10: Infection and vaccination numbers for the Republic of Ireland
starting from January 18th 2021

Our model requires a specific upper bound for the control functions in order
for the optimisation objective to be defined, as seen in Equation 3. One
choice for that upper bound could be the average of the daily vaccinations
over the time period we are studying. However, it is evident from Figures 9(b)
and 10(b) that there is a great fluctuation in the number of vaccines that
were administered which can be attributed to the limited availability of
vaccines at the start of 2021. This leads us to assume that the number of
vaccines administered each day was the maximum of available vaccines for that
day. With that in mind, we designed a time varying upper bound, using each
daily vaccination number as the upper bound for the control function for that
day. Additionally, while the data provides the total number of vaccines to be
administered to the population each day, our model accounts for two age
groups, each of them requiring a separate upper bound. To address this issue,
we apply complimentary percentages of the daily vaccination numbers as upper
bounds to each control function. Specifically, and in agreement with the
policy taken by the Irish government, we consider a higher percentage of the
vaccines to be made a available to the over 65 group. We choose 80% of the
available vaccinations to be used for the over-65 group and let $V_{No}(t)$ be
the number of vaccines administered at time (day) $t$. Then, for time $t$, the
control functions are bounded as follows:

$u_{O}(t)S_{O}(t)<0.8\cdot V_{No}(t)\text{ and }u_{Y}(t)S_{Y}(t)<0.2\cdot
V_{No}(t)$

where each control function is multiplied by the corresponding susceptible
population.

We next set the mean number of people infected by an infectious person, namely
the R0 number. Our model includes 4 age-specific R0 numbers associated with
the $\beta$ parameters in our model, according to Equation (4). In reality,
the R0 number greatly fluctuates depending on individual behaviours and
government restrictions. For that reason, we explore three different choices
for the set of age-specific R0 numbers. These choices are listed in Table 3
along with references to the respective plots.

$R0_{OO}$ | $R0_{YY}$ | $R0_{OY}$ | $R0_{YO}$ | Average $R0$ | Plot before peak | Plot after peak
---|---|---|---|---|---|---
1 | 1 | 0.3 | 0.3 | 1.3 | Figure 11 | Figure 12
1 | 1 | 0.9 | 0.9 | 1.9 | Figure 15 | Figure 16
1.2 | 1.2 | 0.3 | 0.3 | 1.5 | Figure 13 | Figure 14
Table 3: Chosen values for the $R0$ numbers (first four columns). The fifth
column corresponds to the weighted average of the age-specific R0 numbers. The
sixth and seventh column refer to the corresponding plots with the comparative
simulation results, using Irish data before and after the peak in the number
of infections respectively.

(a) Number of infections

(b) Number of vaccinations

Figure 11: Comparison of the infection and vaccination numbers for the
Republic of Ireland starting from January 1st 2021\. The R0 numbers used are
equal to 1 within each age group and 0.3 between age groups.

(a) Infected individuals

(b) Vaccines administered

Figure 12: Comparison of the infection and vaccination numbers for the
Republic of Ireland starting from January 18th 2021.The R0 numbers used are
equal to 1 within each age group and 0.3 between age groups.

(a) Number of infections

(b) Number of vaccinations

Figure 13: Comparison of the infection and vaccination numbers for the
Republic of Ireland starting from January 1st 2021\. The R0 numbers used are
equal to 1.2 within each age group and 0.3 between age groups.

(a) Infected individuals

(b) Vaccines administered

Figure 14: Comparison of the infection and vaccination numbers for the
Republic of Ireland starting from January 18th 2021.The R0 numbers used are
equal to 1.2 within each age group and 0.3 between age groups.

(a) Number of infections

(b) Number of vaccinations

Figure 15: Comparison of the infection and vaccination numbers for the
Republic of Ireland starting from January 1st 2021\. The R0 numbers used are
equal to 1 within each age group and 0.9 between age groups.

(a) Infected individuals

(b) Vaccines administered

Figure 16: Comparison of the infection and vaccination numbers for the
Republic of Ireland starting from January 18th 2021.The R0 numbers used are
equal to 1 within each age group and 0.9 between age groups.

From Figures 11 to 14 it is evident that the optimal vaccination strategy
obtained from our model is successful in quickly reducing the number of
infections even in the case of an average R0 number of $1.9$. This value for
R0 is unrealistic, especially for January 2021, when Ireland was under strict
lockdown measures. However, in this extreme case, while the total number of
infections is high compared to the real numbers, as expected, the decline of
the infection curve is steep, as can be seen in Figure 16. This implies that
even if the country was not under total lockdown during these days, the
optimal vaccination strategy would be able to bring the infections under
control within a short amount of time.

The successful control of infections for all choices of the $R0$ numbers can
be attributed to the distribution of the vaccines to both age groups. The
policy adopted by the Irish government mainly entailed administering vaccines
to the over 65 population as are the higher risk group. However, it is evident
from our simulation that a policy where a portion of the vaccines are
additionally made available to the under 65 group can be more beneficial.

It is interesting to note that in every simulation, the total number of daily
vaccinations resulting from the optimisation are often less than the real
vaccination numbers, i.e. the daily maximum number of available vaccines. This
is a very encouraging observation, considering the implication that less
vaccines can be used to achieve comparable or even better results. It is a
particularly significant result if the optimisation is targeted towards
countries that have gained limited access to vaccines.

##  6 Discussion

Since the end of February 2020, the Republic of Ireland has been affected by
the COVID-19 pandemic, like most countries in the world. For the greater part
of the year no effective treatment or vaccine was present, leaving the
government with one tool to attempt to flatten the curve of infection: the
application of various levels of restrictions, depending on the occurrence of
new daily cases. In light of the vaccines that became available near the end
of 2020, a need arose for efficient vaccination plans based on the
availability of the vaccines and the risk levels of different groups of the
population.

In this work, we introduced a novel compartmental model which more adequately
describes the dynamics of the virus compared to the approaches taken thus far
in the literature. We added compartments representing the stages of infection
and vaccination, and assumed two different age groups sharing the same model
structure. Using this model as a baseline to describe the evolution of the
pandemic in the Republic of Ireland, we applied Optimal Control methods to
obtain a suggested vaccination strategy that minimises the number of
infections under certain restrictions such as the available vaccines. We
simulated the evolution of the pandemic with and without the vaccination
strategy in place for two different scenarios, indicative of strict and loose
restrictions respectively. The main conclusions we drew were:

  1. 1.

The optimal control strategy in both scenarios successfully flattens the curve
of infectious individuals, ensuring that less people will get infected by the
virus due to protection from the vaccine and, most importantly, that less
people will be simultaneously infected thus avoiding exhausting the national
healthcare system.

  2. 2.

In the case of low transmission rates, the vaccination not only flattens the
curve but also significantly reduces the total number of infections,
specifically from 80% to 11.3% for the over 65 group and from 79% to 5.38% for
the under 65 group. However, the pandemic lasts a long time, meaning that a
big part of the population remains in the susceptible compartment, not yet
vaccinated but also not infected due to the low transmission rates, as opposed
to being moved to the protected or recovered compartments.

  3. 3.

When the transmission rates are high however, the vaccination reduces the
number of infections by a much smaller percentage, namely from 99.995% to
85.29% for the over 65 group and from 99.993% to 99.1% for the under 65 group.
However it still succeeds in flattening the curve and ensuring less
simultaneous infections, hence preventing a healthcare crisis. The pandemic in
this case ends in a very short period of time with the biggest part of the
population having been infected and eventually recovered.

  4. 4.

For both scenarios, the optimal strategy suggests focusing on vaccinating the
older population in higher percentages first, while simultaneously vaccinating
part of the younger population. This is a different approach to the one taken
by many states, including Ireland, where the older population gets exclusively
vaccinated first.

All of the above lead us to conclude that an approach involving a combination
of strict and loose measures would be ideal, while the vaccination programme
is taking place. That would ensure that the infection doesn’t spread to such a
high percentage of the population while the restrictions are not as severe so
as to make the situation unbearable for a certain length of time.

To verify our findings, we compared the results of our simulations to the real
infection numbers that occurred during the beginning of 2021 in Ireland, when
the vaccination process was starting across the country. We used three
different sets of transmission numbers, each of them representative of
different levels of restrictions. Our method was successful in bringing the
infection numbers under control within a short amount of time and with a need
for less resources than what was used. This result is particularly important
as it stresses the importance of correct allocation of the available resources
over the quantity of those resources. This means that even countries with
limited access to vaccines are able to effectively reduce the number of
infections.

While our approach led us to a lot of interesting conclusions, there are
certain drawbacks. Our model assumes that all of the parameters remain
constant for the duration of the simulations which is not a realistic
assumption. For example the transmission rates, which were the source of the
vastly different results in our two simulations, do not remain constant for
long periods of time. Furthermore, our model is a deterministic model which
means that there is no accounting for any uncertainty. However it would be
relatively simple to propagate uncertainty through the model by repeatedly
running it with parameter values sourced from an appropriate probability
distribution.

In order to avoid the drawbacks described above, a stochastic model would have
to be introduced or a model with variable parameters. A study on a time-
varying stochastic SEIR model for the control of the Ebola virus can be found
in lekone2006statistical . In addition, as we mentioned earlier, it is
interesting to explore alternative expressions for the objective functional to
be minimised. For example, different weights to the control functions for each
age group can be introduced, modelling the cost of bringing the vaccine to
different populations, or weight factors relating to the mortality rates of
each group, to the numbers of infectious individuals. Moreover, the
introduction of more compartments could enrich the model further, for instance
a compartment for the hospitalised individuals or one for the deceased.
Finally, notable extensions would be the introduction of more age groups and
high risk groups, as well as multiple vaccines and their effectiveness. We
leave these topics as future research topics which we pursue elsewhere.

## Acknowledgments

This publication has emanated from research supported in part by a research
grant from Science Foundation Ireland (SFI) under Grant Number 16/RC/3872 and
is co-funded under the European Regional Development Fund. This work was also
supported by a Science Foundation Ireland COVID-19 Rapid Response Grant number
20/COV/0081. Andrew Parnell’s work was additionally supported by: a Science
Foundation Ireland Career Development Award (17/CDA/4695); an investigator
award (16/IA/4520); a Marine Research Programme funded by the Irish
Government, co-financed by the European Regional Development Fund (Grant-Aid
Agreement No. PBA/CC/18/01); European Union’s Horizon 2020 research and
innovation programme InnoVar under grant agreement No 818144; SFI Centre for
Research Training in Foundations of Data Science 18CRT/6049, and SFI Research
Centre award 12/RC/2289_P2.

##  Appendix A Table of notations

Model states | Symbol
---|---
Susceptible not yet vaccinated (o65 and y65) | $S_{O},S_{Y}$
Received vaccine, waiting for it to take effect (o65 and y65) | $V_{O},V_{Y}$
Received vaccine but was not effective (o65 and y65) | $N_{O},N_{Y}$
Susceptible, refusing or unable to receive vaccine (o65 and y65) | $U_{O},U_{Y}$
Exposed (o65 and y65) | $E_{O},E_{Y}$
Infectious (o65 and y65) | $I_{O},I_{Y}$
Recovered (o65 and y65) | $R_{O},R_{Y}$
Protected from vaccine (o65 and y65) | $P_{O},P_{Y}$
Total number of people in age group (o65 and y65) | $T_{O},T_{Y}$
Model Parameters | Symbol
Rate at which an o65 person infects an o65 person | $\beta_{OO}$
Average number of o65 people infected by an o65 person | $R0_{OO}$
Rate at which a y65 person infects an o65 person | $\beta_{YO}$
Average number of o65 people infected by a y65 person | $R0_{YO}$
Rate at which an o65 person infects a y65 person | $\beta_{OY}$
Average number of y65 people infected by an o65 person | $R0_{OY}$
Rate at which a y65 person infects a y65 person | $\beta_{YY}$
Average number of y65 people infected by a y65 person | $R0_{YY}$
Rate at which exposed becomes infected | $\gamma_{E}$
Rate at which infected becomes recovered | $\gamma_{I}$
Rate at which vaccinated becomes protected | $\gamma_{V}$
Vaccine effectiveness | $\alpha_{V}$
Percentage of over 65s refusing the vaccine | $r_{O}$
Percentage of under 65s refusing the vaccine | $r_{Y}$
Control functions | Symbol
Percentage of over 65s to get vaccinated at time $t$ | $u_{O}(t)$
Percentage of under 65s to get vaccinated at time $t$ | $u_{Y}(t)$
Optimisation Functions and Parameters | Symbol
Hamiltonian function | H
Cost function | $\mathcal{F}$
Age specific weight constants | $W_{O},W_{Y}$

##  Appendix B Existence of an optimal solution

The state equations (2) can be written in compact form as

$g(t,x,u):=\frac{dx(t)}{dt}$ (7)

where $x$ denotes the variable vector
$x=\left(\right.S_{O},V_{O},N_{O},U_{O},E_{O},I_{O},R_{O},P_{O},S_{Y},V_{Y},N_{Y},U_{Y},$
$E_{Y},I_{Y},R_{Y},P_{Y}\left.\right)$. Firstly, we assume $g$ satisfies the
following:

  1. (a)

$|g(t,x,u)|\leq C_{1}(1+|x|+|u|)$

  2. (b)

$|g(t,x,u)-g(t,x^{\prime},u)|\leq C_{2}|x-x^{\prime}|(1+|u|)$

for some positive constants $C_{1},C_{2}$ and for all $t,x,x^{\prime}$ and
$u\in\Omega$. As stated in Fleming1975 , when $g$ is of the class $C^{1}$,
meaning that it is differentiable with continuous first order partial
derivatives, then the above conditions are satisfied by applying suitable
bounds on the partial derivatives of $g$. It is possible to show that such
bounds can be applied in the case of our epidemic model, from the boundedness
of the control and state variables ensuring that Assumptions (a)-(b) are true.

Consider a performance index of the form:

$\mathcal{F}(u(t))=\int_{0}^{T}f(t,x(t),u(t))dt$ (8)

The following theorem, provides the sufficient conditions for the existence of
a solution to the optimal control problem.

###### Theorem B.1

Fleming1975  Suppose that assumptions a-b hold, that $f$ is continuous and
moreover that:

  1. 1.

The set of solutions to the system equations (2) with corresponding control
functions in $\Omega$ is nonempty.

  2. 2.

$\Omega$ is convex.

  3. 3.

$f$ is convex on $\Omega$ and $g$ can be written as $g(t,x,u)=a(t,x)+b(t,x)u$.

Then there exists an optimal control $u^{*}(t)$ with corresponding optimal
trajectory $x^{*}(t)$ minimizing $\mathcal{F}(u(t))$.

We now prove the conditions of the theorem are satisfied for our model. As
noted earlier, conditions (a)-(b) are satisfied. The performance index
integrand is

$f(t,x(t),u(t))=I_{O}(t)+I_{Y}(t)+\frac{W_{O}}{2}u_{O}^{2}(t)+\frac{W_{Y}}{2}u_{Y}^{2}(t)$
(9)

which is continuous, and $\Omega=\left\\{u(t)\in L^{2}(O,T)^{2}\|a\leq
u_{O}(t),u_{Y}(t)\leq b,t\in[0,T]\right\\}$. We now examine conditions 1-3.

  1. 1.

From the Picard-Lindelöf theorem coddington1955theory , if the system
equations are continuous and Lipschitz and the solutions to the system
equations are bounded, then there is a unique solution corresponding to every
admissible control in $\Omega$. It is trivial to verify the continuity of the
state equations (2). Additionally, every element of the variable vector $x(t)$
is bounded by the corresponding age-group totals $T_{O},T_{Y}$. In other
words, $x(t)\in[O,\max\\{T_{O},T_{Y}\\}]^{16}$ which is a compact and convex
set. Finally, the Lipschitz property can be verified with the use of the
following lemma from birkhoff1978ordinary .

###### Lemma 1

If the vector function $\mathbf{X}(\mathbf{x},t)$ is of class $C^{1}$ in a
compact convex domain D, then it satisfies the Lipschitz condition there.

  2. 2.

From the Heine-Borel theorem, it follows that $\Omega$ is convex because as a
closed and bounded set in $\mathbb{R}^{2}$ rudin1976principles .

  3. 3.

The state equations are linearly dependent on the control functions $u_{O}(t)$
and $u_{Y}(t)$ so they can be written in the form $g(t,x,u)=a(t,x)+b(t,x)u$.
The convexity of the integrand $f$ in the objective functional, follows from
the fact that it is quadratic in the control functions.

## References

  * (1) Comirnaty.  European Medicines Agency, https://www.ema.europa.eu/en/medicines/human/EPAR/comirnaty, 2020.  Accessed: 2021-03-16.
  * (2) Coronavirus: Man treated in Dublin hospital as officials trace contacts.  The Irish Times, https://www.irishtimes.com/news/health/coronavirus-man-treated-in-dublin-hospital-as-officials-trace-contacts-1.4189306, 2020\.  Accessed: 2021-03-12.
  * (3) Coronavirus outbreak: Who raises risk warning to very high.  The Irish Times, https://www.irishtimes.com/news/health/coronavirus-outbreak-who-raises-risk-warning-to-very-high-1.4188313, 2020\.  Accessed: 2021-03-12.
  * (4) COVID-19 Vaccine AstraZeneca.  European Medicines Agency, https://www.ema.europa.eu/en/medicines/human/EPAR/covid-19-vaccine-astrazeneca, 2020\.  Accessed: 2021-03-16.
  * (5) COVID-19 Vaccine Janssen.  European Medicines Agency, https://www.ema.europa.eu/en/medicines/human/EPAR/covid-19-vaccine-janssen, 2020\.  Accessed: 2021-04-13.
  * (6) COVID-19 Vaccine Moderna.  European Medicines Agency, https://www.ema.europa.eu/en/medicines/human/EPAR/covid-19-vaccine-moderna#overview-section, 2020\.  Accessed: 2021-03-16.
  * (7) Ireland’s COVID-19 Data Hub.  Government of Ireland, https://covid19ireland-geohive.hub.arcgis.com/pages/detailed-profile-of-cases, 2020\.  Accessed: 2021-03-22.
  * (8) Just 6 % of people will refuse COVID-19 vaccine, says Ipsos MRBI tracker survey for IPHA.  Irish Pharmaceutical Healthcare Association, https://www.ipha.ie/just-6-of-people-will-refuse-covid-19-vaccine-says-ipsos-mrbi-tracker-survey-for-ipha/, 2020\.  Accessed: 2021-04-15.
  * (9) Naming the coronavirus disease (COVID-19) and the virus that causes it.  WHO, https://www.who.int/emergencies/diseases/novel-coronavirus-2019/technical-guidance/naming-the-coronavirus-disease-(covid-2019)-and-the-virus-that-causes-it, 2020\.  Accessed: 2021-03-11.
  * (10) Novel coronavirus - China.  WHO, https://www.who.int/csr/don/12-january-2020-novel-coronavirus-china/en/, 2020\.  Accessed: 2021-03-11.
  * (11) Pfizer and BioNTech announce vaccine candidate against COVID-19 achieved success in first interim analysis from phase 3 study.  Pfizer, https://www.pfizer.com/news/press-release/press-release-detail/pfizer-and-biontech-announce-vaccine-candidate-against, 2020\.  Accessed: 2021-03-16.
  * (12) Pneumonia of unknown cause - China.  WHO, https://www.who.int/csr/don/05-january-2020-pneumonia-of-unkown-cause-china/en/, 2020\.  Accessed: 2021-03-11.
  * (13) Manuel Adrian Acuña-Zegarra, Saúl Díaz-Infante, David Baca-Carrasco, and Daniel Olmos Liceaga.  Covid-19 optimal vaccination policies: A modeling study on efficacy, natural and vaccine-induced immunity responses.  Mathematical Biosciences, page 108614, 2021.
  * (14) Roy M Anderson and Robert M May.  Population biology of infectious diseases: Part I.  Nature, 280(5721):361–367, 1979.
  * (15) Nicolas Bacaër.  A short history of mathematical population dynamics.  Springer Science & Business Media, 2011.
  * (16) Lindsey R. Baden, Hana M. El Sahly, Brandon Essink, et al.  Efficacy and safety of the mrna-1273 sars-cov-2 vaccine.  New England Journal of Medicine, 384(5):403–416, 2020.  PMID: 33378609.
  * (17) Norman TJ Bailey et al.  The mathematical theory of infectious diseases and its applications.  Charles Griffin & Company Ltd, 5a Crendon Street, High Wycombe, Bucks HP13 6LE., 1975.
  * (18) Horst Behncke.  Optimal control of deterministic epidemics.  Optimal Control Applications and Methods, 21(6):269–285, 2000.
  * (19) Garrett Birkhoff and Gian-Carlo Rota.  Ordinary differential equations.  John Wiley & Sons, 1978.
  * (20) V.G. Boltyanski, R.V. Gamkrelidze, E.F. Mishchenko, and L.S. Pontryagin.  The maximum principle in the theory of optimal processes of control.  IFAC Proceedings Volumes, 1(1):464 – 469, 1960.  1st International IFAC Congress on Automatic and Remote Control, Moscow, USSR, 1960.
  * (21) Christopher S Bowman, Julien Arino, and Seyed M Moghadas.  Evaluation of vaccination strategies during pandemic outbreaks.  Mathematical Biosciences & Engineering, 8(1):113, 2011.
  * (22) John Charles Butcher and Nicolette Goodwin.  Numerical methods for ordinary differential equations, volume 2.  Wiley Online Library, 2008.
  * (23) D Cereda, M Tirani, F Rovida, et al.  The early phase of the COVID-19 outbreak in Lombardy, Italy, 2020\.
  * (24) Jasper Fuk-Woo Chan, Shuofeng Yuan, et al.  A familial cluster of pneumonia associated with the 2019 novel coronavirus indicating person-to-person transmission: a study of a family cluster.  The Lancet, 395(10223):514–523, 2021/03/11 2020.
  * (25) Gerardo Chowell, Stefano M Bertozzi, M Arantxa Colchero, et al.  Severe respiratory disease concurrent with the circulation of h1n1 influenza.  New England journal of medicine, 361(7):674–679, 2009.
  * (26) Earl A Coddington and Norman Levinson.  Theory of ordinary differential equations.  Tata McGraw-Hill Education, 1955.
  * (27) Yaqing Fang, Yiting Nie, and Marshare Penny.  Transmission dynamics of the COVID-19 outbreak and effectiveness of government interventions: A data-driven analysis.  Journal of Medical Virology, 92(6):645–659, 2020.
  * (28) Neil M. Ferguson, Derek A. T. Cummings, Christophe Fraser, et al.  Strategies for mitigating an influenza pandemic.  Nature, 442(7101):448–452, 2006.
  * (29) Wendell Fleming and Raymond Rishel.  Deterministic and Stochastic Optimal Control.  Springer-Verlag New York, 1975.
  * (30) Holly Gaff and Elsa Schaefer.  Optimal control applied to vaccination and treatment strategies for various epidemiological models.  Mathematical biosciences & engineering, 6(3):469, 2009.
  * (31) Salisu M. Garba, Jean M.-S. Lubuma, and Berge Tsanou.  Modeling the transmission dynamics of the covid-19 pandemic in south africa.  Mathematical Biosciences, 328:108441, 2020.
  * (32) Alberto Godio, Francesca Pace, and Andrea Vergnano.  SEIR modeling of the Italian epidemic of SARS-CoV-2 using computational swarm intelligence.  International Journal of Environmental Research and Public Health, 17(10):3535, 2020.
  * (33) David S. Hui, Esam I Azhar, et al.  The continuing 2019-nCoV epidemic threat of novel coronaviruses to global health ; The latest 2019 novel coronavirus outbreak in Wuhan, China.  International Journal of Infectious Diseases, 91:264–266, 2021/03/11 2020.
  * (34) Maged N. Kamel Boulos and Estella M. Geraghty.  Geographical tracking and mapping of coronavirus disease covid-19/severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) epidemic and associated events around the world: how 21st century GIS technologies are supporting the global fight against outbreaks and epidemics.  International Journal of Health Geographics, 19(1):8, 2020.
  * (35) Donald E Kirk.  Optimal control theory: an introduction.  Courier Corporation, 2004.
  * (36) Denise Kirschner, Suzanne Lenhart, and Steve Serbin.  Optimal control of the chemotherapy of HIV.  Journal of mathematical biology, 35(7):775–792, 1997.
  * (37) Adam J Kucharski, Timothy W Russell, Charlie Diamond, et al.  Early dynamics of transmission and control of COVID-19: a mathematical modelling study.  The lancet infectious diseases, 20(5):553–558, 2020.
  * (38) Ernest Bruce Lee and Lawrence Markus.  Foundations of optimal control theory.  Robert E. Kreiger Publishing Co., Malabar, Florida, 1967.
  * (39) Sunmi Lee, Michael Golinski, and Gerardo Chowell.  Modeling Optimal Age-Specific Vaccination Strategies Against Pandemic Influenza.  Bulletin of Mathematical Biology, 2012.
  * (40) Phenyo E Lekone and Bärbel F Finkenstädt.  Statistical inference in a stochastic epidemic seir model with control intervention: Ebola as a case study.  Biometrics, 62(4):1170–1177, 2006.
  * (41) Suzanne Lenhart and John T Workman.  Optimal control applied to biological models.  CRC press, 2007.
  * (42) Qun Li, Xuhua Guan, Peng Wu, et al.  Early transmission dynamics in wuhan, china, of novel coronavirus–infected pneumonia.  New England Journal of Medicine, 382(13):1199–1207, 2020.  PMID: 31995857.
  * (43) Qun Li, Xuhua Guan, Peng Wu, et al.  Early Transmission Dynamics in Wuhan, China, of Novel Coronavirus–Infected Pneumonia.  New England Journal of Medicine, 382(13):1199–1207, 2020.  PMID: 31995857.
  * (44) Tao Liu, Jianxiong Hu, Min Kang, et al.  Transmission dynamics of 2019 novel coronavirus (2019-ncov).  bioRxiv, 2020.
  * (45) Hongzhou Lu, Charles W. Stratton, and Yi-Wei Tang.  Outbreak of pneumonia of unknown etiology in wuhan, china: The mystery and the miracle.  Journal of Medical Virology, 92(4):401–402, 2020.
  * (46) Damjan Manevski, Nina Ružić Gorenjec, Nataša Kejžar, and Rok Blagus.  Modeling covid-19 pandemic using bayesian analysis with application to slovene data.  Mathematical Biosciences, 329:108466, 2020.
  * (47) A. G. M’kendrick.  Applications of Mathematics to Medical Problems.  Proceedings of the Edinburgh Mathematical Society, 1925.
  * (48) Amber L Mueller, Maeve S McNamara, and David A Sinclair.  Why does COVID-19 disproportionately affect older people?  Aging (Albany NY), 12(10):9959, 2020.
  * (49) Lennon Ó Náraigh and Áine Byrne.  Piecewise-constant optimal control strategies for controlling the outbreak of covid-19 in the irish population.  Mathematical Biosciences, 330:108496, 2020.
  * (50) Lucien W Neustadt.  Optimization: a theory of necessary conditions.  Princeton University Press, 1976.
  * (51) Calistus N. Ngonghala, Enahoro Iboi, Steffen Eikenberry, Matthew Scotch, Chandini Raina MacIntyre, Matthew H. Bonds, and Abba B. Gumel.  Mathematical assessment of the impact of non-pharmaceutical interventions on curtailing the 2019 novel coronavirus.  Mathematical Biosciences, 325:108364, 2020.
  * (52) Katsuhiko Ogata.  Modern control engineering.  Prentice hall, 2010.
  * (53) T.V. Padma.  Covid vaccines to reach poorest countries in 2023—despite recent pledges.  Nature, 595(7867):342–343, 2021.
  * (54) Rajan Patel, Ira M Longini Jr, and M Elizabeth Halloran.  Finding optimal vaccination strategies for pandemic influenza using genetic algorithms.  Journal of theoretical biology, 234(2):201–212, 2005.
  * (55) Liangrong Peng, Wuyue Yang, Dongyan Zhang, Changjing Zhuge, and Liu Hong.  Epidemic analysis of COVID-19 in China by dynamical modeling.  medRxiv, 2020.
  * (56) Fernando P. Polack, Stephen J. Thomas, Nicholas Kitchin, et al.  Safety and Efficacy of the BNT162b2 mRNA Covid-19 Vaccine.  New England Journal of Medicine, 383(27):2603–2615, 2020.  PMID: 33301246.
  * (57) Lev Semenovich Pontryagin.  Mathematical theory of optimal processes.  Routledge, 2018.
  * (58) Isabelle J. Rao and Margaret L. Brandeau.  Optimal allocation of limited vaccine to control an infectious disease: Simple analytical conditions.  Mathematical Biosciences, 337:108621, 2021.
  * (59) Walter Rudin et al.  Principles of mathematical analysis, volume 3.  McGraw-hill New York, 1976.
  * (60) Oluwaseun Sharomi and Tufail Malik.  Optimal control in epidemiology.  Annals of Operations Research, 251(1-2):55–71, 2017.
  * (61) Mingwang Shen, Zhihang Peng, Yanni Xiao, and Lei Zhang.  Modelling the epidemic trend of the 2019 novel coronavirus outbreak in China.  bioRxiv, 2020.
  * (62) Merryn Voysey, Sue Ann Costa Clemens, Shabir A Madhi, et al.  Safety and efficacy of the ChAdox1 nCov-19 vaccine (AZD1222) against SARS-CoV-2: an interim analysis of four randomised controlled trials in Brazil, South Africa, and the UK.  The Lancet, 397(10269):99–111, 2020.
  * (63) Fan Wu, Su Zhao, Bin Yu, et al.  A new coronavirus associated with human respiratory disease in china.  Nature, 579(7798):265–269, 2020.
  * (64) Joseph T Wu, Kathy Leung, and Gabriel M Leung.  Nowcasting and forecasting the potential domestic and international spread of the 2019-ncov outbreak originating in wuhan, china: a modelling study.  The Lancet, 395(10225):689–697, 2020.
  * (65) Shi Zhao, Qianyin Lin, Jinjun Ran, et al.  Preliminary estimation of the basic reproduction number of novel coronavirus (2019-ncov) in china, from 2019 to 2020: A data-driven analysis in the early phase of the outbreak.  International journal of infectious diseases, 92:214–217, 2020\.
  * (66) Na Zhu, Dingyu Zhang, et al.  A Novel Coronavirus from Patients with Pneumonia in China, 2019\.  New England Journal of Medicine, 382(8):727–733, 2020.  PMID: 31978945.