# Optimal Kernel Tuning Parameter Prediction using Deep Sequence Models

Khawir Mahmood, Jehandad Khan, Hammad Afzal K. Mahmood, H. Afzal are with the
National University of Sciences and Technology, Islamabad, Pakistan.J. Khan is
with AMD Inc., Austin, Texas, USA

###### Abstract

GPU kernels have come to the forefront of computing due to their utility in
varied fields, from high-performance computing to machine learning. A typical
GPU compute kernel is invoked millions, if not billions of times in a typical
application, which makes their performance highly critical. Due to the unknown
nature of the optimization surface, an exhaustive search is required to
discover the global optimum, which is infeasible due to the possible
exponential number of parameter combinations. In this work, we propose a
methodology that uses deep sequence-to-sequence models to predict the optimal
tuning parameters governing compute kernels. This work considers the
prediction of kernel parameters as a sequence to the sequence translation
problem, borrowing models from the Natural Language Processing (NLP) domain.
Parameters describing the input, output and weight tensors are considered as
the _input language_ to the model that emits the corresponding kernel
parameters. In essence, the model _translates_ the problem parameter language
to kernel parameter language. The core contributions of this work are: a)
Proposing that a sequence to sequence model can accurately learn the
performance dynamics of a GPU compute kernel b) A novel network architecture
which predicts the kernel tuning parameters for GPU kernels, c) A constrained
beam search which incorporates the physical limits of the GPU hardware as well
as other expert knowledge reducing the search space. The proposed algorithm
can achieve more than 90% accuracy on various convolutional kernels in MIOpen,
the AMD machine learning primitives library. As a result, the proposed
technique can reduce the development time and compute resources required to
tune unseen input configurations, resulting in shorter development cycles,
reduced development costs, and better user experience.

###### Index Terms:

Convolutional Neural Networks, Deep Learning, GPU Compute Kernels, Neural
Machine Translation, Sequence Models, Tuning Parameters

##  I Introduction

General-Purpose Computing on Graphics Processing Units (GPGPU) has enabled
modern Machine Learning (ML) algorithms to perform at-par and, in some cases,
superior to humans [1, 2, 3]. ML libraries such as Tensorflow [4] and Pytorch
[5] convert the data flow graph into smaller operations, which are optimized
to exploit the parallelism offered by a GPU. These operations take the form of
compute-intensive routines called _kernels_ ; performance optimized kernels
require a deep understanding of the target hardware to extract the maximum
possible performance and enable an algorithm to run at high efficiency. The
design process creates a tight coupling between the kernel and the target
hardware, making it challenging to innovate and adapt to changing
requirements.

The standard approach to enhance code porting and re-usability between
different variants of the same hardware class is to parameterize the kernels
and perform a parameter search on the target hardware to ascertain their
optimum values to reach peak performance. In some cases, these kernel
parameters are hand-tuned for a particular architecture by engineers with a
deep understanding of the hardware architecture, resulting in a tedious
procedure that is difficult to scale for production software.

For production software, search techniques discover output kernel tuning
parameters wherein kernels with different parameter values are compiled and
benchmarked to arrive at the optimal parameter values [6, 7]. Various
heuristics are employed to avoid performing an exhaustive search, but the cost
in terms of computing resources and time remains significant. This approach
makes it difficult to generalize to unseen input configurations, a critical
feature for production software.

Predicting the optimum tuning parameters is a challenging task due to the non-
trivial nature of operations involved (such as convolutions and batch
normalization), hardware complexity, and code design. These factors manifest
as a non-linear optimization surface that is discrete and inexpressible as an
analytical model; furthermore, the combinatorial nature of the parameters
makes it infeasible to traverse the search space exhaustively. This paper
proposes an ML-based approach using sequence models (Recurrent Neural
Networks) for predicting these parameters. By incorporating a novel constraint
satisfaction procedure, the proposed method ensures that the predicted
parameters are always valid for the target architecture as well as improve
performance by reducing the search space. Moreover, the Recurrent Neural
Network (RNN) learns the underlying dynamics of the kernel and can generalize
to unseen problem configurations.

The core contributions of this work are: a) Proposing that a sequence to
sequence model can accurately learn the performance dynamics of a GPU compute
kernel b) A novel network architecture which predicts the kernel tuning
parameters for GPU kernels, c) A constrained beam search which incorporates
the physical limits of the GPU hardware as well as other expert knowledge
reducing the search space.

The paper’s organization is as follows: Section II reviews existing kernel
parameter tuning techniques, while Section IV details the mechanics of the
proposed algorithm. Section V further improves the Seq2Seq framework by
incorporating a constrained beam search. Section VI compares the accuracy of
the proposed work against classical machine learning methods as well as
contrasts the performance of the various proposed architectures, section VII
concludes the paper.

##  II Prior Work

Prior art for kernel tuning consists of two main areas: model-based, and
black-box approaches. Model-based approaches aim to develop an analytical
model to inform the search process about the behavior of the target hardware.
On the other hand, model-less methods tend to treat the problem as an
optimization problem aimed to discover the optimum by using domain-specific
variants of some known optimization algorithms. The following paragraphs
highlight the salient works in these areas and contrast them with the proposed
approach. While these approaches succeed to varying degrees in designing
efficient heuristics to discover the optimum tuning parameters, the actual
optimum may only be found by an exhaustive evaluation of the search space, as
is typically the case with combinatorial optimization problems.

The lack of an underlying model in black-box optimization-based approaches
enables these techniques on various hardware targets without modification.
However, this also means that the algorithm needs to be run from the start for
each program, implying results from one execution may not improve successive
runs. On the other hand, model-based approaches suffer from the limitations of
the models they rely on, which may be different for various architectures.
Machine-Learning based approaches promise to strike a balance between the two
by automating the model design process as well as the ability to generalize
based on the supplied data. In the following paragraphs, we begin by reviewing
prevalent work in model-free approaches to kernel tuning, followed by a few
methods based on model-based techniques.

Papenhausen et al. [8] extend the ant colony optimization technique [9] to
discover optimal tuning parameters for kernels generated as part of the PPCG
compiler framework [10]. By incorporating performance into the cost model and
sharding the search space, the algorithm can efficiently search the parameter
space. Workhoven et al. [6] present a framework that enables the tuning of
various kernels with different optimization algorithms. Kernel tuner
contributes implementations for various optimization algorithms such as Basin
Hopping, Differential Evolution, and the Firefly algorithm [11]. CLTune by
Nugetern et al. [7] and AutoTuner by Bruel et al. [12] present frameworks
that, while do not introduce new techniques, enable tuning for kernels using
search strategies such as simulated annealing and particle swarm optimization.

In contrast to optimization-based approaches, model-based methods aim to
develop analytical models to discover the ideal set of parameters for a given
kernel and hardware combination. Analytical models have the advantage of being
faster to evaluate (compared to searching a space); however, they may lead to
some losses due to model approximations and mismatches. Lim et al. [13] use
static analysis as the driving force to collect information about a GPU
program, augmenting it with dynamic data collected during program execution.
An analytical model, after that, predicts the set of optimum parameters based
on the collected data. _RoofLine_ models [14] predict the performance limits
of a kernel and identify whether a kernel is compute-bound or memory-bound.
However, translating such information to concrete performance parameter
recommendations might not always be possible. The model makes certain
assumptions about program execution, which may or may not be accurate in all
cases. For example, in specific compute-intensive kernels, it might be better
to launch fewer wave-fronts (or warps) and have each thread perform an
increased amount of work. Furthermore, each new device generation obsoletes
older models, making them less ideal for production use.

Machine Learning (ML) methods discover the underlying model using execution
data, which promises to strike a balance between searching and hand-tuned
models. Bergstra et al. [15] employ a boosted regression tree to predict the
execution time of a compute kernel; this allows a search algorithm to discard
bad candidates and search in a smaller space quickly. Falch et al. [16] adopt
a similar approach; they use an Artificial Neural Network to directly prune
the search space by randomly sampling the parameter space to train a neural
network. SnowPack [17] and StaTuner [18] use support vector machines as the
regression models to estimate the runtime of a kernel and avoid benchmarking
kernels by launching them.

The approaches discussed above use the machine learning model as an estimator
of the kernel runtime; in contrast, DeepTune [19] adopts a method that is
closest to work presented here. In DeepTune, a deep neural network learns
features directly from the kernel code like Natural Language Processing (NLP)
models, using them to predict the kernel performance parameters. The proposed
work improves upon this work by removing the need to parse the source code and
incorporate hardware constraints dictated by the hardware and the kernel
design. Furthermore, the proposed algorithm can predict as many as eight
performance parameters as compared to only one in DeepTune.

##  III Problem Description

Modern neural network models rely on various operations to form compute graphs
which learn the target function, however, the most useful one, particularly
for vision-based tasks is the convolution operation. The convolution operation
accounts for a considerable chunk of the training time spent to learn its
parameters. Therefore, improved performance parameters for convolutional
layers would yield considerable dividends in terms of better performance. Due
to their criticality in modern networks and tuning complexity, they are the
ideal candidates for a performance parameter search study.

###  III-A The convolution operation

Convolution in machine learning refers to the traditional de-convolution
operation in engineering, where the learning algorithm learns the weights of
the kernel using the data it is presented [20]. The important notion here is
that learning happens over the weights of the convolutional kernel as opposed
to traditional machine vision where the weights of the kernel are known
_apriori_. Therefore, the convolution operation takes two inputs; the image
being convolved and the weights of the kernels and produces the convolved
image as the output. Formally:

$\displaystyle O$ $\displaystyle=I\circledast W$ (1)

Where $O$ represents the output image and $I$ and $W$ represent the input
image and the convolution weights, respectively. The convolution operation
itself is represented as $\circledast$. Keep in mind, the notion of input and
output changes in during back-propagation with $W$ and $I$ being outputs
during the backward weights and backward data operations. However, both those
can again be expressed as regular convolution and therefore are irrelevant to
the present discussion.

To improve computational efficiency and enhance throughput multiple input
images are grouped together in _batches_ and convolved together. Therefore,
$I$ and $W$ become four dimensional tensors instead of being two dimensional
images with various channels (i.e., three dimensional tensors). Therefore,
equation 1 may be written as

$\displaystyle O(n,k,h_{o},w_{o})$
$\displaystyle=I(n,c,h_{i},w_{i})\circledast W(k,c,y,x)$ (2)

Where $O$, $I$ and $W$ are the output, input, and weight tensors respectively
as before. The letters in parenthesis represent the dimensionality of each
tensor with $n$ representing the batch-size, $h_{o}$ and $w_{o}$ represent the
height and width of the output image, while $h_{i}$ and $w_{i}$ is the height
and width of the input image. $k$ represents the number of channels in the
output image, $c$ represents the number of channels in the input image, $y$
and $x$ is the height and width of the convolution filter, respectively.

The dimensions of a tensor collectively are known as the _Tensor Descriptor_
for a particular tensor. Therefore, the input tensor is described by the
tensor descriptor containing $n,c,h_{i},w_{i}$, collectively these parameters
represent the feature space for our kernel performance parameter predictor. A
computer kernel (whether on the GPU or a CPU) converts the input and weight
hyper-volume to an output hyper-volume, with the amount of memory and number
of computation required directly proportional to their collective hyper-
volume.

A kernel developer employs various optimizations around memory access
patterns, data layouts, memory hierarchy, access width etc. to compute the
output volume as efficiently as possible. Collectively, these choices comprise
the design space of the kernel author, with different points in this space
resulting in different performance characteristics. The difficulty lies in
choosing the ideal parameters for a given set of input tensor descriptors,
since the relative volume of the input tensors deeply affects the choice of
performance parameters to achieve optimum performance. The proposed algorithm
predicts these algorithmic choices based on the various input tensor
descriptors. The choice of convolution as the problem space is motivated by
the complexity and importance of this operation in various machine learning
scenarios.

##  IV Method

The proposed model is based on traditional language models that translate a
sequence of token in one language to an equivalent sequence of words in
another language. These models are called Sequence-to-Sequence models [21].
One class of such models is called Encoder-Decoder models due to the division
of the network architecture in two distinct networks connected by a hidden
state [22]. For the task at hand, such models perform remarkably better
compared to traditional classification networks in the following ways: a) A
single model can learn to predict all the tuning parameters of a kernel, thus
learning the joint probability distribution over the target domain. b) Due to
the recurrent nature of the networks, they are better able to learn and
generalize the correlations between different parameters. However, there are
certain disadvantages to this choice of architecture as well; for instance,
due to the token-based nature of the input and the output, the network cannot
learn to predict unseen ranges.

The integral nature of the input, output, and weight tensor descriptors makes
it tempting to use their numerical values to predict the integral tuning
parameters. However, there are certain caveats associated with this thinking.
If an estimator function learns to predict based on the integral value of the
height of an image, the learned function would be a continuous one, making
fractional image sizes valid – which is not correct. Furthermore, the behavior
of the kernel may change at specific values, such as when some hardware
resource is exhausted, which again is difficult to learn in a continuous
setting. Therefore, in practice, the authors observed that encoding the tensor
descriptors as well as the predicted performance parameters as one-hot encoded
tokens yielded better results, as evidenced by the results section. In
contrast to typical language processing tasks that embed the input sequence of
tokens in some space [23], [24], this work performs an end-to-end conversion
between input tokens representing characteristics of the convolution operation
and tokens representing performance parameters.

###  IV-A The Basic Encoder-Decoder Model

In their seminal work Sutskever, et al. [21] used the RNN architecture to
outperform state-of-the-art on sequence-to-sequence prediction tasks. The
strength of the model lies in its ability to train a single end-to-end model
on the source and target sequences and to handle variable-length input and
output sequences. The encoders and decoders comprise of Long Short-Term
Memory(LSTM) [21] or Gated Recurrent Unit (GRU) [25] cells.

Figure 1 shows the encoder-decoder model unrolled in the time-dimension; the
encoder and decoder network itself consists of only one cell. For each sample,
the input sequence $[x_{0},x_{1},...x_{t}]$ is fed to the encoder in
consecutive time-steps $0\rightarrow t$. The encoder updates its internal
state $e$ at each time step, i.e., $e_{0}\rightarrow
e_{1}\rightarrow....e_{t}$ using information from the input at that time step
$i$ while preserving the information it has already learned from previous
time-steps $0\rightarrow i-1$. Thus, at final time-step _t_ , the whole input
sequence is represented by the internal state $e_{t}$ also known as thought
vector $w$.

Based on $e_{t}$, at each time step $t$, the decoder generates the output
sequence $[y\hat{}_{0},y\hat{}_{1},...y\hat{}_{t}]$ as a probability
distribution over the output parameter range. Note that, at each time-step
after $t=0$, the predicted value from previous time step $y\hat{}_{t-1}$ is
also used for making the prediction, i.e., updating the decoder’s internal
state $d_{t}$.

Figure 1: Encoder-Decoder Architecture

###  IV-B Attention

The Encoder-Decoder model relies on a single fixed-length _thought vector_ to
represent the complete source sequence while discarding the useful information
in the intermediate encoder states, resulting in inefficient handling of
longer sequences. Augmenting the Encoder-Decoder networks with an _Attention_
mechanism [26, 27] enables the decoder to focus on different sub-sequences
during each step, overcoming the limitations discussed above.

This work proposes a modified version of _Bahdanau-attention_ [26] by not
reusing intermediate outputs at each time-step, likely because of independent
output tokens. The model shown in Figure 2 consists of a bi-directional
encoder LSTM called _pre-attention bi-LSTM_ , an _attention mechanism_ , and a
decoder LSTM called _post-attention LSTM_.

Figure 2: Proposed Attention Model

For each sample, the input sequence $[x_{0},x_{1},...x_{t}]$ is fed to the
pre-attention bi-LSTM in consecutive time-steps $0\rightarrow t$. For each
input $x_{i}$, the pre-attention bi-LSTM outputs the forward
$\overrightarrow{a_{i}}$ and backward $\overleftarrow{a_{i}}$ activations
(internal states), concatenated together form $a_{i}$.

$a_{i}=\\{\overrightarrow{a_{i}},\overleftarrow{a_{i}}\\}$ (3)

At each time-step $t$ the attention weights $\alpha_{(i,t^{\prime})}$ over all
input time-steps $t^{\prime}$ are calculated using the previous post-attention
LSTM internal state $s_{i-1}$ and concatenated activations $a_{i}$ of the pre-
attention bi-LSTM. A Neural Network consisting of a single _Dense_ layer
calculates the weights $d$ which are then passed through a _Softmax_ layer to
get a probability distribution of attention weights $\alpha$.

$d_{(i,t^{\prime})}=\textup{Dense}(s_{i-1},a_{t^{\prime}})\alpha_{(i,t^{\prime})}=\frac{\exp(d_{(i,t^{\prime})})}{\sum_{k=1}^{t}\exp(d_{(i,t^{\prime})})}$
(4)

The attention weights $\alpha_{(i,t^{\prime})}$ are combined together to form
the Context Vector $C_{i}$ which is calculated by taking a sum of the
attention weights $\alpha_{(i,t^{\prime})}$ weighted by the respective input
activations $a_{t^{\prime}}$.

$C_{i}=\sum_{j=1}^{t}\alpha_{(i,t^{\prime})}a_{t^{\prime}}$ (5)

The Context Vector $C_{i}$ is used by the post-attention LSTM to sequentially
predict $y_{i}$ as a probability distribution over the output parameter range.

###  IV-C Adding Convolutions - Hybrid Model

While the attention module, coupled with bi-directional LSTMs capture long
term dependencies. However, Convolutional Neural Networks (CNN) are better
adept at capturing spatial relationships in data and may be computed more
efficiently than RNNs [28, 29]. Popular due to their success in image-based
classification tasks, convolutional layers have been used successfully in many
other applications such as image captioning [30], change detection [31], and
video activity recognition [32, 33, 34].

In this spirit, the authors implemented a hybrid model shown in Figure 3 that
uses convolutional (Conv) layers as the encoder followed by bi-LSTMs for the
decoder.

Figure 3: Conv1D encoder and Bi-LSTM decoder architecture

The complete input sequence passes through multiple _1-dimensional
Convolutional layers_ ($Conv1D$). Each layer consists of various filters $f$,
of size $k$ that perform convolutions with stride $s$. On an input sequence of
length $t$, a $Conv1D$ forms $f$ representations of the data, each of size
$o=\frac{t-k}{s}+1$. Traditionally each convolutional layer is followed by a
pooling layer. However, in the case of kernel parameter prediction, adding
pooling layers reduces model performance, likely due to the reduced size of
the feature space, in contrast to image processing, where images have a large
number of pixels. The output from the last $Conv1D$ layer is flattened to form
a 1-dimensional vector $d=f\times o$, used to predict the complete output
sequence. The decoder consists of multiple Bi-LSTM layers that sequentially
predict the output sequence based on the flattened input from the Conv
Encoder. The Bi-LSTM layer helps map interdependencies in the output sequence.
Finally, a Dense layer is used with Softmax activation to create the output
sequence likelihood with the highest value corresponding to the predicted
output.

Figure 4 shows a variation of the Hybrid Model, where only the last states
(forward $\overrightarrow{b_{t}}$ and reverse $\overleftarrow{b_{0}}$) from
bi-LSTM-1 are passed on to bi-LSTM-2 removing false temporal dependencies
formed amongst predicted output sequences when using multiple bi-LSTM layers.

Figure 4: Conv1D encoder and layered bi-LSTM decoder architecture

##  V Constrained Beam Search

The decoder of an RNN computes the token probability at the time $i$
conditioned on the hidden vector state (post attention) as well as all output
tokens until time $t-1$, implying that an erroneous prediction at any time may
result in an incorrect prediction. This problem may be mitigated using a
technique known as _Beam Search_ and is commonly used to improve prediction
performance in case of sequence to sequence models [35, 36]. Instead of only
propagating the word with the highest probability at each time step, the
method propagates $k$ subsequences with the highest likelihood at each time
step as shown in Figure 5. The parameter $k$ is known as _Beam Width_.
Incorporating beam search in the models presented above resulted in
considerable gains in prediction accuracy, as discussed in detail below.

Figure 6 depicts how beam search fits in the proposed architecture, where
$y_{0\rightarrowtail t}$ denotes time-steps or output parameters, and
$y^{0\rightarrowtail k}$ denotes the top $k$ choices. Thus, at each time-step
$i$, $k$ best choices are kept based on the complete input sequence $x$ and
the values leading up to the current time-step $y_{0\rightarrowtail t-1}$.

Figure 5: Beam Search For selecting $k$ values with highest probabilities at
each step Figure 6: Change in shape of output kernel parameter sequence after
applying Beam Search

###  V-A Applying Output Parameter Constraint Satisfaction

The performance parameters for a compute kernel translate to the number of
used hardware resources, such as the number of registers, amount of shared
memory, types of instructions, etc. These resources have physical limits
which, if violated, may either render a kernel invalid or result in drastic
performance cliffs. This domain knowledge is invaluable in terms of reducing
the search space for the predicted optimal tuning parameter. While the notion
of constraining the output sequence in language translation models is not
useful, for the case at hand presents a unique opportunity to improve
prediction performance. These constraints are in the form of Boolean predicate
functions that take the tensor descriptor and the predicted sub-sequence as
inputs and results in the feasibility of the predicted parameter combination.
Discarding the candidates that violate the output parameter constraints,
considerably improves accuracy. Therefore, at each time-step the $k$ candidate
sub-sequences are sure to be valid for the kernel in question, Figure 7
depicts this process in more detail.

Figure 7: The overall flow of Improved Hybrid _Seq2Seq_ Framework showing
language encoding, prediction, beam-search, constraint-satisfaction, and
language decoding pipeline

##  VI Results and Analysis

Multiple kernels from MIOpen [37] were selected as candidates, to ascertain
the efficacy of the proposed neural architecture. MIOpen is AMD’s open-source
library for high-performance machine learning primitives supporting different
operations such as convolution, batchnorm, and RNNs. MIOpen ships with optimal
performance parameters for these kernels for various convolution
configurations, these input parameters were processed according to the
techniques discussed earlier. Table I outlines the selected kernels as well as
the number and range of their performance parameters; the total number of
possible parameters is vast, making it a formidable task to predict these
parameters accurately.

TABLE I: Kernel Output Parameters and Ranges Kernel | Parameters | Range/Values
---|---|---
ConvAsm1x1U | read_size | 1-4
k_mult | 1,4,8,16,32
chunks_per_wave | 1-16
chunk_size | 1,2,4,8,16,32,64
n_mult | 1-8
c_mult | 1,2,4,8,16,32
waves_c_in_group | 1-8
waves_k_in_group | 1,2,4,8
ConvOclDirectFwd1x1 | grp_tile1 | $2^{0-4}$
grp_tile0 | $2^{0-8}$
in_tile1 | $2^{0-5}$
in_tile_0 | $2^{0-5}$
out_pix_tile1 | 0,1
out_pix_tile0 | 0,1,2,4
n_out_pix_tiles | $2^{0-6}$
n_in_data_tiles | $2^{0-11}$
n_stacks | 0,1
ConvAsmBwdWrW1x1 | read_size | 1-4
c_per_gpr | 1,2,4,8,16
c_mult | 1,2,4,8,16
k_per_gpr | 1,2,4,8,16
k_mult | 1,2,4,8,16
n_per_gpr | 1,2,4,8,16
n_part_cnt | 1-8
chunk_size | 1,2,4,8,16
short_store | 0,1
data_prefetch | 0-4
ConvAsmBwdWrW3x3 | limit_wave_cnt | 0-9
reverse_inout | 0,1
chunk_size | 8,16
k_per_wave | 1,2,4,8
pipe_lines_depth | 1-16
n_per_group | 1-8

To better assess the performance of various techniques, it is useful to define
the following metrics:

The _mean of the prediction accuracy_ of each output parameter is defined as

$AverageAccuracy={\frac{1}{n}}\sum_{i=1}^{n}(\frac{C_{i}}{S_{Test}}*100)$

Where $n$ is the total number of output parameters for the kernel, $C_{i}$ is
the number of correct predictions of parameter $i$ and $S_{Test}$ is the total
number of test samples. The percentage of _perfect predictions_ in the test
sample, where the complete predicted output parameter set matches the optimum
parameter set for the test sample is given by:-

$PerfectPrediction=\frac{count(p_{0\rightarrow t}=a_{0\rightarrow
t})}{S_{Test}}\times 100$

Where $p_{0\rightarrow t}$ is the predicted output parameter/sequence set,
$a_{0\rightarrow t}$ is the actual output parameter/sequence set, and
$S_{Test}$ is the total number of test samples. This metric quantifies a joint
prediction of the parameter sequence. Therefore, it is a more useful metric
for evaluating model performance since the key objective is to predict the
complete set of optimum kernel tuning parameters.

Classical ML techniques such as regression trees, clustering, and multi-layer
perceptrons set a baseline for prediction performance. This exercise provided
insight into the data and assurance that no trivially discoverable patterns
were present. Results from this exercise are presented in Table II.

TABLE II: Average Accuracy kernel parameter prediction benchmarks using Classic ML techniques Kernel | Precision | Samples | D.Tree | GNB | k-NN | SVC | RFC | Average | Best
---|---|---|---|---|---|---|---|---|---
ConvAsm1x1U | Full | 7,038 | 61.52 | 48.45 | 52.18 | 52.16 | 64.80 | 54.60 | 64.80
Half | 3,291 | 78.25 | 58.15 | 62.53 | 61.14 | 80.34 | 66.43 | 80.34
ConvOclDirectFwd1x1 | Full | 21,124 | 77.58 | 68.41 | 68.77 | 69.32 | 79.02 | 71.92 | 79.02
Half | 3,349 | 81.49 | 72.87 | 76.58 | 76.35 | 82.19 | 77.06 | 82.19
ConvAsmBwdWrW1x1 | Full | 5,100 | 71.00 | 66.00 | 70.71 | 71.43 | 74.96 | 70.02 | 74.96
Half | 1,606 | 73.04 | 61.37 | 66.46 | 66.09 | 75.34 | 67.28 | 75.34
ConvAsmBwdWrW3x3 | Full | 698 | 60.63 | 47.10 | 52.66 | 55.07 | 65.70 | 54.71 | 65.70
Half | 447 | 67.04 | 51.11 | 60.37 | 68.15 | 72.22 | 61.67 | 72.22

Table III details the performance metrics of the models described in Section
IV for different kernels and precisions. Sequential models consistently
perform better than the classic ML techniques as measured by the Average
Accuracy metric. The hybrid layered decoder model described in Section IV
outperforms all other variants in both average and perfect predictions.

TABLE III: Average Accuracy and Perfect Prediction Results of Seq2Seq Model Variations Kernel | Precision | Enc-Dec | Attn | Attn-2 | Hybrid | Hybrid-2
---|---|---|---|---|---|---
Avg | Pft | Avg | Pft | Avg | Pft | Avg | Pft | Avg | Pft
ConvAsm1x1U | Full | 55.27 | 4.26 | 55.47 | 2.13 | 56.71 | 1.85 | 63.25 | 6.96 | 62.59 | 6.53
Half | 79.85 | 31.21 | 80.57 | 27.27 | 80.87 | 30.00 | 82.16 | 31.21 | 82.23 | 31.52
ConvOclDirectFwd1x1 | Full | 82.84 | 15.76 | 82.33 | 15.05 | 82.78 | 16.75 | 82.63 | 15.52 | 83.64 | 18.93
Half | 87.06 | 30.15 | 86.07 | 25.37 | 85.57 | 26.27 | 87.10 | 29.25 | 86.73 | 31.34
ConvAsmBwdWrW1x1 | Full | 75.39 | 8.24 | 74.00 | 4.31 | 73.63 | 4.51 | 75.96 | 10.39 | 75.78 | 10.00
Half | 77.08 | 12.53 | 75.52 | 11.18 | 74.72 | 12.42 | 76.58 | 8.07 | 77.76 | 13.66
ConvAsmBwdWrW3x3 | Full | 73.43 | 11.59 | 65.46 | 7.25 | 72.46 | 11.94 | 75.24 | 12.86 | 75.48 | 15.43
Half | 76.67 | 17.78 | 77.41 | 13.33 | 75.93 | 11.11 | 77.78 | 17.78 | 79.26 | 22.22
Avg: Average Accuracy, Pft: Prefect Prediction, Enc-Dec: Encoder Decoder,
Attn: Attention
Attn-2: Modified Attention, Hybrid-2: Hybrid Layered Decoder Model

Augmenting the Hybrid Seq2Seq Layered Decoder Model with a constrained beam
search ensures that a valid set of output parameter sequences is produced and
improves results prediction metrics, as shown in Table IV. Although it
increases the number of outputs from a single output kernel parameter set to
$k$ sets (where $k$ is the beam width), compared to empirical optimization
techniques, even keeping $k=100$ has practical significance due to an increase
in Perfect Prediction metric.

TABLE IV: Results of Hybrid Seq2Seq Framework with output parameter constraint satisfaction fused Beam Search Kernel | Precision | Metric | Hybrid-2 with Beam Search | Hybrid-2 with OPCS-BS
---|---|---|---|---
k=10 | k=30 | k=50 | k=100 | k=10 | k=30 | k=50 | k=100
ConvAsm1x1U | Full | Avg | 79.79 | 87.48 | 90.04 | 92.97 | 79.99 | 88.23 | 91.16 | 94.07
Pft | 19.60 | 28.55 | 33.24 | 39.20 | 19.89 | 30.40 | 34.52 | 40.91
Half | Avg | 94.58 | 97.58 | 98.48 | 99.02 | 95.27 | 98.30 | 99.20 | 99.62
Pft | 67.88 | 76.36 | 80.30 | 82.73 | 68.48 | 78.79 | 81.52 | 87.27
ConvOclDirectFwd1x1 | Full | Avg | 96.90 | 98.73 | 98.93 | 99.33 | 97.20 | 98.82 | 99.42 | 99.77
Pft | 69.85 | 87.17 | 91.81 | 94.94 | 72.50 | 90.49 | 93.94 | 95.79
Half | Avg | 97.65 | 99.20 | 99.44 | 99.73 | 97.68 | 99.27 | 99.50 | 99.80
Pft | 79.10 | 85.97 | 88.36 | 91.04 | 79.70 | 86.57 | 89.25 | 91.64
ConvAsmBwdWrW1x1 | Full | Avg | 90.04 | 95.39 | 96.65 | 97.92 | 90.12 | 96.22 | 97.84 | 98.71
Pft | 34.90 | 47.45 | 54.31 | 65.10 | 36.47 | 52.94 | 60.00 | 71.37
Half | Avg | 92.05 | 96.83 | 97.33 | 98.63 | 92.17 | 97.39 | 99.19 | 99.50
Pft | 40.37 | 57.76 | 62.73 | 67.08 | 44.10 | 62.11 | 65.84 | 73.29
ConvAsmBwdWrW3x3 | Full | Avg | 94.29 | 98.10 | 99.29 | 99.52 | - | - | - | -
Pft | 58.57 | 81.43 | 87.14 | 91.43 | - | - | - | -
Half | Avg | 96..67 | 99.63 | 100.00 | 100.00 | - | - | - | -
Pft | 77.78 | 88.89 | 91.11 | 95.56 | - | - | - | -
Hybrid-2: Hybrid Layered Decoder Model, OPCS-BS: Output Parameter Constraint
Satisfaction fused Beam Search
Avg: Average Accuracy, Pft: Perfect Prediction

The RNN based model used here has various hyperparameters, which have a
considerable impact on prediction accuracy. Key findings from trying various
combinations for different models presented in Section IV are as follows:

  * •

The size of the internal state vector $|e|$ for the encoder-decoder
architecture strongly influences the quality of results, $|e|=256$ gives the
best results. Moreover, setting decoder LSTM _Dropout_ and _Recurrent Dropout_
to 0.2 (i.e., 20%) improved test results by preventing over-fitting.

  * •

The tunable hyperparameters of both attention and modified attention models
include the number of input/output layers, cell sizes of the pre-attention bi-
LSTM $n_{a}$, post-attention LSTM $n_{s}$, and the number of nodes in the
Dense layer $n_{d}$. $n_{a}=256$, $n_{s}=512$, and $n_{d}=2$ give the best
results.

  * •

In addition to the number of Conv and Sequential layers, tunable
hyperparameters for the Hybrid Model include those specifically related to
Conv1D as well as bi-LSTM. Hyperparameters of Conv1D include the number of
filters $f$, the size of each filter, also known as kernel size $k$, and the
convolutional stride $s$. For the bi-LSTM, hyperparameters include cell size
of each unit, dropout, and recurrent dropout. A cell size of _256_ for most of
the kernels gives the best results.

###  VI-A Training convergence analysis

Analysis of train-test accuracy improvement (and correspondingly loss
reduction) with the increasing number of epochs gives insight into the
learning of the Hybrid Seq2Seq model on various kernels. The change in
accuracy and loss as a function of the number of epochs is shown for all four
Kernels in Figures 8,9,10,11. ConvAsm1x1U and ConvOclDirectFwd1x1 kernels in
half-precision mode display closest training and test accuracy improvement
(and corresponding loss reduction) with successive epochs. This learning
performance corresponds to the comparatively better results seen for these
kernels. Overall, additional data would aid learning for certain full-
precision kernels.

Figure 8: Change in Accuracy and Loss w.r.t Number of Epochs for ConvAsm1x1U
Half and Full Precision Figure 9: Change in Accuracy and Loss w.r.t Number of
Epochs for ConvOclDirectFwd1x1 Half and Full Precision Figure 10: Change in
Accuracy and Loss w.r.t Number of Epochs for ConvAsmBwdWrW1x1 Half and Full
Precision Figure 11: Change in Accuracy and Loss w.r.t Number of Epochs for
ConvAsmBwdWrW3x3 Half and Full Precision

##  VII Conclusion

_Seq2Seq_ models are powerful ML constructs that perform well in many domains.
In this work, Seq2Seq models in the presence of limited data samples have
shown the capacity to learn a discrete and complex search space with more than
95% accuracy for many kernels. Convolutional (Conv) Layers are most suitable
for forming data abstractions and representations, especially for problems
such as kernel parameter prediction where the sequences are not related in the
temporal domain but have dependencies that originate from complex functions
with multiple dimensions. This work indicates that it is possible to learn the
computing behavior of a GPU using an RNN and use this model to predict the
optimal kernel tuning parameters for compute kernels.

## References

  * [1] C. Li, Y. Yang, M. Feng, S. T. Chakradhar, and H. Zhou, “Optimizing memory efficiency for deep convolutional neural networks on GPUs,” _CoRR_ , vol. abs/1610.03618, 2016. [Online]. Available: http://arxiv.org/abs/1610.03618
  * [2] A. V. Sinitskiy and V. S. Pande, “Physical machine learning outperforms “human learning” in Quantum Chemistry,” _arXiv e-prints_ , p. arXiv:1908.00971, Aug 2019.
  * [3] D. Silver, A. Huang, C. J. Maddison, A. Guez, L. Sifre, G. van den Driessche, J. Schrittwieser, I. Antonoglou, V. Panneershelvam, M. Lanctot, S. Dieleman, D. Grewe, J. Nham, N. Kalchbrenner, I. Sutskever, T. Lillicrap, M. Leach, K. Kavukcuoglu, T. Graepel, and D. Hassabis, “Mastering the game of Go with deep neural networks and tree search,” _Nature_ , vol. 529, pp. 484–503, 2016. [Online]. Available: http://www.nature.com/nature/journal/v529/n7587/full/nature16961.html
  * [4] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard _et al._ , “TensorFlow: A system for large-scale machine learning,” in _12th $\\{$USENIX$\\}$ Symposium on Operating Systems Design and Implementation ($\\{$OSDI$\\}$ 16)_, 2016, pp. 265–283.
  * [5] A. Paszke, S. Gross, S. Chintala, and G. Chanan, “PyTorch: Tensors and dynamic neural networks in python with strong GPU acceleration,” 2017, https://pytorch.org/.
  * [6] B. van Werkhoven, “Kernel tuner: A search-optimizing GPU code auto-tuner,” _Future Generation Computer Systems_ , vol. 90, pp. 347 – 358, 2019. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0167739X18313359
  * [7] C. Nugteren and V. Codreanu, “CLTune: A generic auto-tuner for OpenCL kernels,” in _2015 IEEE 9th International Symposium on Embedded Multicore/Many-core Systems-on-Chip_ , Sep. 2015, pp. 195–202.
  * [8] E. Papenhausen and K. Mueller, “Coding ants: Optimization of GPU code using ant colony optimization,” _Computer Languages, Systems & Structures_, vol. 54, pp. 119–138, 05 2018.
  * [9] M. Dorigo, M. Birattari, and T. Stutzle, “Ant colony optimization,” _IEEE computational intelligence magazine_ , vol. 1, no. 4, pp. 28–39, 2006.
  * [10] S. Verdoolaege, J. Carlos Juega, A. Cohen, J. Ignacio Gomez, C. Tenllado, and F. Catthoor, “Polyhedral parallel code generation for CUDA,” _ACM Transactions on Architecture and Code Optimization (TACO)_ , vol. 9, no. 4, pp. 1–23, 2013.
  * [11] X.-S. Yang and X. He, “Firefly algorithm: recent advances and applications,” _arXiv preprint arXiv:1308.3898_ , 2013.
  * [12] P. Bruel, A. Gonzalez, Marcos, and A. Goldman, “Autotuning CUDA compiler parameters for heterogeneous applications using the OpenTuner framework,” _Concurrency and Computation Practice and Experience_ , 03 2017.
  * [13] R. V. Lim, B. Norris, and A. D. Malony, “Autotuning GPU kernels via static and predictive analysis,” _CoRR_ , vol. abs/1701.08547, 2017. [Online]. Available: http://arxiv.org/abs/1701.08547
  * [14] E. Konstantinidis and Y. Cotronis, “A quantitative roofline model for GPU kernel performance estimation using micro-benchmarks and hardware metric profiling,” _Journal of Parallel and Distributed Computing_ , vol. 107, pp. 37–56, 04 2017.
  * [15] J. Bergstra, N. Pinto, and D. Cox, “Machine learning for predictive auto-tuning with boosted regression trees,” in _2012 Innovative Parallel Computing (InPar)_ , May 2012, pp. 1–9.
  * [16] T. L. Falch and A. C. Elster, “Machine learning based auto-tuning for enhanced OpenCL performance portability,” _CoRR_ , vol. abs/1506.00842, 2015. [Online]. Available: http://arxiv.org/abs/1506.00842
  * [17] R. Sing, P. Wood, R. Gupta, S. Bagchi, and I. Laguna, “Snowpack: Efficient parameter choice for GPU kernels via static analysis and statistical prediction,” in _Proceedings of the 8th Workshop on Latest Advances in Scalable Algorithms for Large-Scale Systems_ , ser. ScalA ’17. New York, NY, USA: ACM, 2017, pp. 8:1–8:8. [Online]. Available: http://doi.acm.org/10.1145/3148226.3148235
  * [18] R. Gupta, I. Laguna, D. H. Ahn, T. Gamblin, S. Bagchi, and F. X. Lin, “STATuner: Efficient tuning of CUDA kernels parameters,” in _Supercomputing Conference (SC 2015)_ , 2015.
  * [19] C. Cummins, P. Petoumenos, Z. Wang, and H. Leather, “End-to-end deep learning of optimization heuristics,” in _2017 26th International Conference on Parallel Architectures and Compilation Techniques (PACT)_ , Sep. 2017, pp. 219–232.
  * [20] Y. LeCun, Y. Bengio _et al._ , “Convolutional networks for images, speech, and time series,” _The handbook of brain theory and neural networks_ , vol. 3361, no. 10, p. 1995, 1995.
  * [21] I. Sutskever, O. Vinyals, and Q. V. Le, “Sequence to sequence learning with neural networks,” _CoRR_ , vol. abs/1409.3215, 2014. [Online]. Available: http://arxiv.org/abs/1409.3215
  * [22] L. Fridman, “MIT deep learning basics: Introduction and overview with TensorFlow,” https://medium.com/tensorflow/mit-deep-learning-basics-introduction-and-overview-with-tensorflow-355bcd26baf0, accessed: 2019-08-18.
  * [23] T. Mikolov, I. Sutskever, K. Chen, G. Corrado, and J. Dean, “Distributed representations of words and phrases and their compositionality,” _CoRR_ , vol. abs/1310.4546, 2013. [Online]. Available: http://arxiv.org/abs/1310.4546
  * [24] J. Pennington, R. Socher, and C. D. Manning, “GloVe: Global vectors for word representation,” in _Empirical Methods in Natural Language Processing (EMNLP)_ , 2014, pp. 1532–1543. [Online]. Available: http://www.aclweb.org/anthology/D14-1162
  * [25] K. Cho, B. van Merrienboer, Ç. Gülçehre, F. Bougares, H. Schwenk, and Y. Bengio, “Learning phrase representations using RNN encoder-decoder for statistical machine translation,” _CoRR_ , vol. abs/1406.1078, 2014. [Online]. Available: http://arxiv.org/abs/1406.1078
  * [26] B. Dzmitry, K. Cho, and Y. Bengio, “Neural machine translation by jointly learning to align and translate,” 2014, cite arxiv:1409.0473Comment: Accepted at ICLR 2015 as oral presentation. [Online]. Available: http://arxiv.org/abs/1409.0473
  * [27] M. Luong, H. Pham, and C. D. Manning, “Effective approaches to attention-based neural machine translation,” _CoRR_ , vol. abs/1508.04025, 2015. [Online]. Available: http://arxiv.org/abs/1508.04025
  * [28] N. Kalchbrenner and P. Blunsom, “Recurrent convolutional neural networks for discourse compositionality,” _CoRR_ , vol. abs/1306.3584, 2013. [Online]. Available: http://arxiv.org/abs/1306.3584
  * [29] X. Zhang, J. Zhao, and Y. LeCun, “Character-level convolutional networks for text classification,” in _Advances in Neural Information Processing Systems 28_ , C. Cortes, N. D. Lawrence, D. D. Lee, M. Sugiyama, and R. Garnett, Eds. Curran Associates, Inc., 2015, pp. 649–657. [Online]. Available: http://papers.nips.cc/paper/5782-character-level-convolutional-networks-for-text-classification.pdf
  * [30] Q. You, H. Jin, Z. Wang, C. Fang, and J. Luo, “Image captioning with semantic attention,” in _The IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2016.
  * [31] W. Luo, W. Liu, and S. Gao, “A revisit of sparse coding based anomaly detection in stacked rnn framework,” in _Proceedings of the IEEE International Conference on Computer Vision_ , 2017, pp. 341–349.
  * [32] W. Zhang, J. Hu, Y. Feng, and Q. Liu, “Information-propogation-enhanced neural machine translation by relation model,” _CoRR_ , 09 2017.
  * [33] M. Elbayad, L. Besacier, and J. Verbeek, “Pervasive attention: 2D convolutional neural networks for sequence-to-sequence prediction,” _CoRR_ , vol. abs/1808.03867, 2018. [Online]. Available: http://arxiv.org/abs/1808.03867
  * [34] J. Gehring, M. Auli, D. Grangier, D. Yarats, and Y. N. Dauphin, “Convolutional sequence to sequence learning,” in _Proceedings of the 34th International Conference on Machine Learning - Volume 70_ , ser. ICML’17. JMLR.org, 2017, pp. 1243–1252. [Online]. Available: http://dl.acm.org/citation.cfm?id=3305381.3305510
  * [35] M. Freitag and Y. AlOnaizan, “Beam search strategies for neural machine translation,” _CoRR_ , vol. abs/1702.01806, 2017. [Online]. Available: http://arxiv.org/abs/1702.01806
  * [36] S. Wiseman and A. M. Rush, “Sequence-to-sequence learning as beam-search optimization,” _CoRR_ , vol. abs/1606.02960, 2016. [Online]. Available: http://arxiv.org/abs/1606.02960
  * [37] J. Khan, P. Fultz, A. Tamazov, D. Lowell, C. Liu, M. Melesse, M. Nandhimandalam, K. Nasyrov, I. Perminov, T. Shah _et al._ , “Miopen: An open source library for deep learning primitives,” _arXiv preprint arXiv:1910.00078_ , 2019.


# Social vs. individual age-dependent costs of imperfect vaccination

Fabio A. C. C. Chalub111Corresponding author<EMAIL_ADDRESS>Center for
Mathematics and Applications (NOVA Math) and Department of Mathematics, NOVA
FCT, Caparica, 2829-516, Portugal  Paulo Doutor  Center for Mathematics and
Applications (NOVA Math) and Department of Mathematics, NOVA FCT, Caparica,
2829-516, Portugal  Paula Patrício  Center for Mathematics and Applications
(NOVA Math) and Department of Mathematics, NOVA FCT, Caparica, 2829-516,
Portugal  Maria do Céu Soares  Center for Mathematics and Applications (NOVA
Math) and Department of Mathematics, NOVA FCT, Caparica, 2829-516, Portugal

###### Abstract

In diseases with long-term immunity, vaccination is known to increase the
average age at infection as a result of the decrease in the pathogen
circulation. This implies that a vaccination campaign can have negative
effects when a disease is more costly (financial or health-related costs) for
higher ages. This work considers an age-structured population transmission
model with imperfect vaccination. Our aim is to compare the social and
individual costs of vaccination, assuming that disease costs are age-
dependent. A model coupling pathogen deterministic dynamics for a population
consisting of juveniles and adults, both assumed to be rational agents, is
introduced. The parameter region for which vaccination has a positive social
impact is fully characterized and the Nash equilibrium of the vaccination game
is obtained. Finally, collective strategies designed to promote voluntary
vaccination, without compromising social welfare, are discussed.

##  1 Introduction

In a voluntary vaccination scheme, in which the vaccine is perceived – truly
or falsely – as risky, herd immunity will never be attained in a population
composed of rational individuals [2]. Just before vaccine coverage reaches the
herd immunity threshold, rational individuals will stop to be vaccinated as
the perceived risk of the vaccine will equal the perceived risk of the
disease, which will be small at this point. Therefore, herd immunity will be
obtained through vaccination only if there are incentives to be vaccinated
(and to vaccinate the dependents) or punishment of non-vaccinated individuals
(e.g., the exclusion of the school system). Since the seminal work [2], other
models considered the coupling between the deterministic disease dynamics with
game-theoretical models for individual decisions within the population, cf.
[11, 19, 3, 8]. See also [5] for a previous work of the present group of
authors in models for voluntary vaccinations in seasonal diseases.

A pathogen in a partially vaccinated population (i.e., below herd immunity
level) will circulate slower than in a non-vaccinated population. Assuming
long-term immunity for recovered individuals, a partial vaccination will be
the increase in the average age of infected individuals [13].

Furthermore, it is naive to expect, for any particular vaccine, a 100%
efficacy, cf. [15]. Depending on the precise details of the disease dynamics
and its effect on the population, an imperfect vaccination scheme may have
adverse collective effects.

Let us see a particular example. Consider a disease in which the effect is
different in juveniles and adults as for chickenpox, rubella, or Zika. The
infection has an overall mild effect in juveniles, but when the virus infects
adults, particularly pregnant women, the health consequences can be more
severe [14, 4, 1]. While full coverage of a perfect vaccine would prevent the
disease spread and a free circulation of a highly infectious virus will
asymptotically turn it into a child disease, with mild economic and health
effects, a partial vaccination may be pernicious. As a consequence, it is
important to find the parameter region where it is better to vaccinate than to
not vaccinate, and also it is important to establish if it is possible to move
continuously in the parameter region such that full coverage can be attained
within acceptable social costs.

Models using game-theoretical arguments for the study of imperfect vaccination
were presented in [10] and [12]. In both cases, three Nash equilibria were
found in the model and the vaccination coverage for the Nash equilibrium may
be higher than for the social optimum, depending on the costs of vaccination.
In the former case, the authors determine whether the optimal vaccination
coverage may be achieved through individual action, comparing two different
vaccination scenarios for chickenpox (USA and Israel). In the latter a model
with reinfection is considered, and two of the three Nash equilibria are
evolutionarily stable, with a catastrophe from the high-vaccination to the
low-vaccination scenario, where the effect of vaccination is worse for the
population as a whole. We will introduce a precise definition of Nash
equilibria in the context of vaccination games shortly; for now, it is enough
to consider a situation in which all individuals in the population
simultaneously and freely minimize the joint cost of the disease and the
vaccine.

In this work, we compare social _vs_. individual interests regarding
vaccination and disease costs and investigate if it is possible to promote
voluntary vaccination and still satisfy both interests. For that, we consider
an age-structured model with age-dependent costs, permanent immunity, and
imperfect vaccination and use a game theory approach to analyze individual
decisions.

We finish the introduction with the outline of the paper. In Section 2, we
introduce the model and present some basic results, including the explicit
expression for the basic reproduction number, and the characterization of
equilibria and their stability. In the sequel, we discuss the model, analyzing
first the social costs of vaccination and then, using techniques from game
theory, the effects of considering voluntary vaccination and individual
interests; in particular, we define Nash equilibrium within the context of the
present work. In Section 3, we present numerical simulations based on typical
values for child diseases to study socially cost-efficient parameters regions,
Nash equilibria of the vaccination games, parameters region such that rational
individuals accept or do not accept to be vaccinated, and how shared costs
between individuals and the society can dramatically influence the endemic
equilibria of the model. We conclude in Section 4 with a summary.

##  2 Methods

###  2.1 The model

We consider an age-structured population divided into two groups: juveniles
and adults. Each individual is vaccinated at birth with probability
$p\in[0,1]$. The vaccine is imperfect, with efficacy $\lambda\in[0,1]$,
meaning that with probability $\lambda$ it confers life-long immunity, while
with probability $1-\lambda$ the immunity only lasts during the juvenile phase
($1/\nu$ yrs). The model diagram is represented in Fig. 1. The relevant set of
values is presented at Table 1, while model variables are defined at Table 2.

Figure 1: Schematic diagram of the SIR model for juveniles and adults. The transition rate between both age groups is given by $\nu$. Vaccination (with coverage $p$) provides long term immunity for a fraction $\lambda$ of the individuals and temporary (i.e., during the juvenile phase) for a fraction $1-\lambda$. Disease transmission $\beta$ and recovering $\gamma$ are assumed to be independent of the age group. Parameter | Description | Value | Unity
---|---|---|---
$\mu>0$ | birth/mortality rate | 1/70 | yrs-1
$\gamma>0$ | recovering rate | 365/12 | yrs-1
$\beta>0$ | transmission rate | such that ${\mathcal{R}}_{0}=8$ | yrs-1 per capita
$\nu>0$ | rate of immunity loss | 1/15 | yrs-1
$p\in[0,1]$ | vaccine coverage |  | non-dimensional
$\lambda\in[0,1]$ | vaccine efficacy |  | non-dimensional
Table 1: Values used in this work. Parameters $\mu$, $\nu$, $\gamma$, and $\beta$ are not disease-specific and were chosen as an illustration in the range of Chickenpox and Rubella that served as motivation [6]. The value of $\beta$ was obtained from Eq. (11) at demographic equilibrium. In Fig. 4 we consider a range of values ${\mathcal{R}}_{0}$. Variable |  Description
---|---
$V$ |  Fraction of individuals vaccinated at birth
$S_{{\mathrm{J}}}$ |  Fraction of susceptible juveniles
$I_{{\mathrm{J}}}$ |  Fraction of infectious juveniles,
$R_{{\mathrm{J}}}$ |  Fraction of juveniles with life-long immunity (due to recovery or vaccination)
$S_{{\mathrm{A}}}$ |  Fraction of susceptible adults
$I_{{\mathrm{A}}}$ |  Fraction of infectious adults
$R_{{\mathrm{A}}}$ |  Fraction of adults with life-long immunity (due to recovery or vaccination)
$N_{{\mathrm{J}}}$ |  Fraction of juveniles (equal to $V+S_{{\mathrm{J}}}+I_{{\mathrm{J}}}+R_{{\mathrm{J}}}$)
$N_{{\mathrm{A}}}$ |  Fraction of adults (equal to $S_{{\mathrm{A}}}+I_{{\mathrm{A}}}+R_{{\mathrm{A}}}$)
Table 2: Compartment variables used in the model; c.f. Eqs. (1)–(7).

The model can be represented by the following system of differential
equations:

$\displaystyle V^{\prime}$ $\displaystyle=\mu p(1-\lambda)N_{{\mathrm{A}}}-\nu
V\ ,$ (1) $\displaystyle S_{{\mathrm{J}}}^{\prime}$
$\displaystyle=\mu(1-p)N_{A}-\nu
S_{{\mathrm{J}}}-\beta(I_{{\mathrm{J}}}+I_{{\mathrm{A}}})S_{{\mathrm{J}}}\ ,$
(2) $\displaystyle I_{{\mathrm{J}}}^{\prime}$
$\displaystyle=\beta(I_{{\mathrm{J}}}+I_{{\mathrm{A}}})S_{{\mathrm{J}}}-\nu
I_{{\mathrm{J}}}-\gamma I_{{\mathrm{J}}}\ ,$ (3) $\displaystyle
R_{{\mathrm{J}}}^{\prime}$ $\displaystyle=\mu p\lambda N_{{\mathrm{A}}}+\gamma
I_{{\mathrm{J}}}-\nu R_{{\mathrm{J}}}\ ,$ (4) $\displaystyle
S_{{\mathrm{A}}}^{\prime}$ $\displaystyle=\nu(V+S_{{\mathrm{J}}})-\mu
S_{{\mathrm{A}}}-\beta(I_{{\mathrm{J}}}+I_{{\mathrm{A}}})S_{{\mathrm{A}}}\ ,$
(5) $\displaystyle I_{{\mathrm{A}}}^{\prime}$
$\displaystyle=\beta(I_{{\mathrm{J}}}+I_{{\mathrm{A}}})S_{{\mathrm{A}}}+\nu
I_{{\mathrm{J}}}-\mu I_{{\mathrm{A}}}-\gamma I_{{\mathrm{A}}}\ ,$ (6)
$\displaystyle R_{{\mathrm{A}}}^{\prime}$ $\displaystyle=\nu
R_{{\mathrm{J}}}+\gamma I_{{\mathrm{A}}}-\mu R_{{\mathrm{A}}}\ .$ (7)

The total population
$N=V+S_{{\mathrm{J}}}+I_{{\mathrm{J}}}+R_{{\mathrm{J}}}+S_{{\mathrm{A}}}+I_{{\mathrm{A}}}+R_{{\mathrm{A}}}$
is constant, and, therefore, we set $N(t)=1$ for all $t\geq 0$. Furthermore,
we define the juvenile and adult population by
$N_{{\mathrm{J}}}:=V+S_{{\mathrm{J}}}+I_{{\mathrm{J}}}+R_{{\mathrm{J}}}$ and
$N_{{\mathrm{A}}}:=S_{{\mathrm{A}}}+I_{{\mathrm{A}}}+R_{{\mathrm{A}}}=1-N_{{\mathrm{J}}}$,
respectively. Adding Eqs. (5), (6) and (7), we conclude that
$N_{{\mathrm{A}}}^{\prime}=\nu(1-N_{{\mathrm{A}}})-\mu N_{{\mathrm{A}}}$. We
say that a population is in _demographic equilibrium_ if $N_{{\mathrm{J}}}$
and $N_{{\mathrm{A}}}$ are constants. In that case

$N_{{\mathrm{J}}}(t)=N_{{\mathrm{J}}}^{*}:=\frac{\mu}{\mu+\nu}\ ,\quad
N_{{\mathrm{A}}}(t)=N_{{\mathrm{A}}}^{*}:=\frac{\nu}{\nu+\mu}\ .$ (8)

Both the disease-free and endemic equilibrium can be readily obtained. Their
stability depends on the value of the critical parameter $\mathcal{R}_{p}$,
obtained using the next generation matrix approach [17]. More explicitly, we
state:

######  Theorem 1.

For any value of $p\in[0,1]$, there is one equilibrium solution of Eqs.
(1)–(7), called _the disease-free solution_ , given by

$\displaystyle V^{{\mathrm{df}}}$
$\displaystyle:=N_{{\mathrm{J}}}^{*}p(1-\lambda),$ $\displaystyle
S_{{\mathrm{J}}}^{\mathrm{df}}:=N_{{\mathrm{J}}}^{*}(1-p),$ $\displaystyle
R_{{\mathrm{J}}}^{\mathrm{df}}$ $\displaystyle:=N_{{\mathrm{J}}}^{*}p\lambda,$
$\displaystyle
S_{{\mathrm{A}}}^{\mathrm{df}}:=N_{{\mathrm{A}}}^{*}(1-p\lambda),$
$\displaystyle R_{{\mathrm{A}}}^{\mathrm{df}}$
$\displaystyle:=N_{{\mathrm{A}}}^{*}p\lambda,$ $\displaystyle
I_{{\mathrm{J}}}^{\mathrm{df}}:=I_{{\mathrm{A}}}^{\mathrm{df}}=0.$

Let the effective reproduction number be

$\displaystyle\mathcal{R}_{p}$
$\displaystyle:=\\!\frac{\beta}{\gamma+\mu}\\!\left[\frac{\mu+\nu+\gamma}{\nu+\gamma}S_{J}^{\mathrm{df}}+S_{A}^{\mathrm{df}}\right]$
(9)
$\displaystyle=\frac{\beta}{\gamma+\mu}\frac{\mu(\mu+\gamma+\nu)(1-p)+\nu(\nu+\gamma)(1-\lambda
p)}{(\gamma+\nu)(\mu+\nu)}\ .$

Then

  * •

If $\mathcal{R}_{p}<1$ the only equilibrium solution of Eqs. (1)–(7) is the
disease-free solution, which is locally asymptotically stable.

  * •

If $\mathcal{R}_{p}>1$ the disease-free solution is unstable. Furthermore,
there is a second equilibrium solution of Eqs. (1)–(7), called _the endemic
solution_ , given by

$\displaystyle V^{\mathrm{en}}$
$\displaystyle:=N_{{\mathrm{J}}}^{*}p(1-\lambda)=\frac{\mu
p(1-\lambda)}{\mu+\nu},$ $\displaystyle S^{\mathrm{en}}_{{\mathrm{J}}}$
$\displaystyle:=\frac{N_{J}^{*}(1-p)\nu}{\nu+\beta
I^{\mathrm{en}}}=\frac{\mu\nu(1-p)}{(\mu+\nu)(\nu+\beta I^{\mathrm{en}})},$
$\displaystyle R^{\mathrm{en}}_{J}$
$\displaystyle:=\frac{\gamma}{\nu}I_{{\mathrm{J}}}^{\mathrm{en}}+N_{{\mathrm{J}}}^{*}p\lambda=N_{J}^{*}\left[\frac{(1-p)\gamma\beta
I^{\mathrm{en}}}{(\gamma+\nu)(\nu+\beta I^{\mathrm{en}})}+p\nu\right],$
$\displaystyle S^{\mathrm{en}}_{{\mathrm{A}}}$ $\displaystyle:=\mu
N_{{\mathrm{A}}}^{*}\frac{(1-p)\nu+p(1-\lambda)(\nu+\beta
I^{\mathrm{en}})}{(\mu+\beta I^{\mathrm{en}})(\nu+\beta I^{\mathrm{en}})}\ ,$
$\displaystyle R_{{\mathrm{A}}}^{\mathrm{en}}$
$\displaystyle:=\frac{\gamma}{\mu}I^{\mathrm{en}}+N_{{\mathrm{A}}}^{*}p\lambda,$
$\displaystyle I^{\mathrm{en}}_{{\mathrm{J}}}$
$\displaystyle:=\frac{N^{*}_{{\mathrm{J}}}(1-p)\beta
I^{\mathrm{en}}\nu}{(\nu+\gamma)(\nu+\beta I^{\mathrm{en}})},$ $\displaystyle
I^{\mathrm{en}}_{{\mathrm{A}}}$ $\displaystyle:=\frac{\mu
N^{*}_{{\mathrm{A}}}\beta
I^{\mathrm{en}}}{\mu+\gamma}\left[\frac{p(1-\lambda)}{\mu+\beta
I^{\mathrm{en}}}+\frac{(1-p)\nu}{(\mu+\beta I^{\mathrm{en}})(\nu+\beta
I^{\mathrm{en}})}\right.$ $\displaystyle\hskip
216.81pt+\left.\frac{(1-p)\nu}{(\nu+\gamma)(\nu+\beta
I^{\mathrm{en}})}\right].$

Finally, the total number of infectious individuals at the endemic equilibrium
is given by

$I^{\mathrm{en}}:=I^{\mathrm{en}}_{{\mathrm{J}}}+I^{\mathrm{en}}_{{\mathrm{A}}}=\frac{b_{1}+\sqrt{b_{1}^{2}+4b_{2}b_{0}}}{2b_{2}}\
,$ (10)

where

$\displaystyle b_{0}$
$\displaystyle:=\mu\nu\left[\beta(\mu(\mu+\gamma+\nu)(1-p)+\nu(\nu+\gamma)(1-\lambda
p))\right.$ $\displaystyle\hskip
216.81pt\left.-(\gamma+\mu)(\gamma+\nu)(\mu+\nu)\right]$
$\displaystyle(b_{0}>0\Leftrightarrow\mathcal{R}_{p}>1)\ ,$ $\displaystyle
b_{1}$ $\displaystyle:=\beta^{2}\nu\mu((\gamma+\nu)(1-\lambda
p)+\mu(1-p))-\beta(\gamma+\mu)(\gamma+\nu)(\mu+\nu)^{2},$ $\displaystyle
b_{2}$ $\displaystyle:=\beta^{2}(\gamma+\mu)(\gamma+\nu)(\mu+\nu)\ .$

###### Proof.

The disease-free solution is immediate after imposing
$I^{\mathrm{df}}_{{\mathrm{J}}}=I^{\mathrm{df}}_{{\mathrm{A}}}=0$ in the
stationary (i.e., ${}^{\prime}=0$) solution of the System (1)–(7). Following
[17], we consider the compartments corresponding to infectious individuals to
be $x=(I_{{\mathrm{J}}},I_{{\mathrm{A}}})$ and the remaining compartments
corresponding to non-infectious classes
$y=(V,S_{{\mathrm{J}}},R_{{\mathrm{J}}},S_{{\mathrm{A}}},R_{{\mathrm{A}}})$.
We define the rate of appearance of new infections as
$\mathcal{F}(x,y)=(\beta(I_{{\mathrm{J}}}+I_{{{\mathrm{A}}}})S_{{\mathrm{J}}},\beta(I_{{\mathrm{J}}}+I_{{{\mathrm{A}}}})S_{{{\mathrm{A}}}}))$
and the remaining transition terms as $\mathcal{V}(x,y)=(\nu I_{J}+\gamma
I_{J},-\nu I_{J}+(\gamma+\mu)I_{{\mathrm{A}}})$. Hence, System (1) can be
written as

$x^{\prime}=\mathcal{F}(x,y)-\mathcal{V}(x,y),\quad y^{\prime}=g(x,y)\ ,$

for an appropriate function $g$. We define the matrices

$F=\left[\frac{\partial\mathcal{F}_{i}}{\partial
x_{j}}(x_{0},y_{0})\right]=\begin{bmatrix}\beta
S^{{\mathrm{df}}}_{{\mathrm{J}}}&\beta S^{{\mathrm{df}}}_{{\mathrm{J}}}\\\
\beta S^{{\mathrm{df}}}_{{\mathrm{A}}}&\beta
S^{{\mathrm{df}}}_{{\mathrm{A}}}\end{bmatrix}$

and

$V=\left[\frac{\partial\mathcal{V}_{i}}{\partial
x_{j}}(x_{0},y_{0})\right]=\begin{bmatrix}\nu+\gamma&0\\\
-\nu&\mu+\gamma\end{bmatrix}\ ,$

where $(x_{0},y_{0})$ represents the disease free equilibrium. It’s
straightforward to verify conditions $(A_{1})$ to $(A_{5})$ of Theorem 2 in
[17], hence we conclude that the effective reproduction number
$\mathcal{R}_{p}$ is given by the spectral radius of the next generation
matrix

$FV^{-1}=\frac{\beta}{(\gamma+\mu)(\gamma+\nu)}\begin{bmatrix}(\gamma+\mu+\nu)S^{{\mathrm{df}}}_{{\mathrm{J}}}&(\gamma+\nu)S^{{\mathrm{df}}}_{{\mathrm{J}}}\\\
(\gamma+\mu+\nu)S^{{\mathrm{df}}}_{\mathcal{A}}&(\gamma+\nu)S^{{\mathrm{df}}}_{\mathcal{A}}\end{bmatrix}\
,$

i.e.,

$\mathcal{R}_{p}:=\frac{\beta}{\gamma+\mu}\left[\frac{\mu+\nu+\gamma}{\nu+\gamma}S_{J}^{{\mathrm{df}}}+S_{A}^{{\mathrm{df}}}\right]\
.$

The stability follows from [17], namely the disease-free equilibrium is
locally asymptotically stable if ${\mathcal{R}_{p}<1}$, and unstable if
$\mathcal{R}_{p}>1$. For the computation of the endemic equilibrium we follow
the same techniques as before; in this case, however, the stationary solution
implicitly depends on the value of $I^{\mathrm{en}}$, the solution of
$\wp(I)=0$, where $\wp(I):=-b_{2}I^{2}+b_{1}I+b_{0}$. ∎

As an immediate consequence of Theorem 1, we write

$\mathcal{R}_{p}=\mathcal{R}_{0}\left[1-p\left(1-\frac{(1-\lambda)\nu(\nu+\gamma)}{\mu(\mu+\gamma+\nu)+\nu(\nu+\gamma)}\right)\right],$

with

${\mathcal{R}}_{0}:=\left.{\mathcal{R}}_{p}\right|_{p=0}=\frac{\beta}{\gamma+\mu}\left[\left(1+\frac{\mu}{\nu+\gamma}\right)N_{J}^{*}+N_{A}^{*}\right]\
.$ (11)

Furthermore,

######  Theorem 2.

Let
$\Gamma=\\{(V,S_{{\mathrm{J}}},I_{{\mathrm{J}}},R_{{\mathrm{J}}},S_{{\mathrm{A}}},I_{{\mathrm{A}}},R_{{\mathrm{A}}}):S_{{\mathrm{J}}}\leq
S_{{\mathrm{J}}}^{\mathrm{df}},S_{{\mathrm{A}}}\leq
S_{{\mathrm{A}}}^{\mathrm{df}},\linebreak V\leq
V^{\mathrm{df}},N_{{\mathrm{A}}}\leq N_{{\mathrm{A}}}^{*}\\}$, and consider
the model given by Eqs. (1)–(7). Then

  * •

If $\mathcal{R}_{p}<1$ the only equilibrium solution of the System (1)–(7) is
the disease-free solution, which is globally asymptotically stable in
$\Gamma$.

  * •

If $\mathcal{R}_{p}>1$ the disease-free solution is unstable. The System
(1)–(7) is uniformly persistence.

###### Proof.

The set $\Gamma$ is positively invariant. Following the notation from the
proof of Thm. 1 we define $f(x,y)=(F-V)x-\mathcal{F}+\mathcal{V}$. We have
that $f(x,y)\geq 0$ with $f(x,y_{0})=0$ in $\Gamma$, $F\geq 0$, $V^{-1}\geq 0$
and $V^{-1}F$ is irreducible. Moreover,
$(0,y)=(0,N^{*}_{{\mathrm{J}}},0,N^{*}_{{\mathrm{A}}},0)$ is a globally
asymptotically stable (GAS) equilibrium of the system $y^{\prime}=g(0,y)$.
Hence, by [16, Thm. 2.2], we conclude that the disease-free solution is GAS in
$\Gamma$ for $\mathcal{R}_{p}<1$ and that, for $\mathcal{R}_{p}>1$, the system
is uniformly persistent. ∎

Finally, it is straightforward to prove that

######  Proposition 3.

Consider

$\lambda>\lambda_{c}:=1-\frac{(\gamma+\mu)(\mu+\nu)}{\beta\nu}$ (12)

and $\mathcal{R}_{0}>1$. Then, there is a critical vaccination coverage

$p_{\mathrm{c}}:=\frac{\mu(\mu+\gamma+\nu)+\nu(\nu+\gamma)}{\mu(\mu+\gamma+\nu)+\lambda\nu(\gamma+\nu)}\left(1-\frac{1}{\mathcal{R}_{0}}\right)\in(0,1)$
(13)

such that for any $p>p_{\mathrm{c}}$ the disease free solution is globally
asymptotically stable in $\Gamma$.

###  2.2 Social cost

At the endemic equilibrium, we define a social cost function (per unit of
time) depending on the disease incidence and disease cost for both juveniles
and adults and on the vaccination costs:

$\displaystyle\phi(p,\lambda)$
$\displaystyle:=c_{{\mathrm{A}}}^{\mathrm{d}}\beta(I^{\mathrm{en}}_{{\mathrm{J}}}+I^{\mathrm{en}}_{{\mathrm{A}}})S^{\mathrm{en}}_{{\mathrm{A}}}+c_{{\mathrm{J}}}^{\mathrm{d}}(\beta(I^{\mathrm{en}}_{{\mathrm{J}}}+I^{\mathrm{en}}_{{\mathrm{A}}})S^{\mathrm{en}}_{{\mathrm{J}}}+\nu
I^{\mathrm{en}}_{{\mathrm{J}}})+c^{\mathrm{v}}\delta\mu pN_{{\mathrm{A}}}^{*}$
$\displaystyle=c_{{\mathrm{A}}}^{\mathrm{d}}(\gamma+\mu)I_{{\mathrm{A}}}^{{\mathrm{en}}}+c_{{\mathrm{J}}}^{\mathrm{d}}(\gamma+\nu)I_{{\mathrm{J}}}^{{\mathrm{en}}}+c^{\mathrm{v}}\delta\mu
pN_{{\mathrm{A}}}^{*}$
$\displaystyle=c_{{\mathrm{A}}}^{\mathrm{d}}[(\gamma+\mu)I_{{\mathrm{A}}}^{{\mathrm{en}}}+\varepsilon(\gamma+\nu)I_{{\mathrm{J}}}^{{\mathrm{en}}}+r\delta\mu
pN_{{\mathrm{A}}}^{*}],$

where $c_{{\mathrm{A}}}^{\mathrm{d}}>0$ and $c_{{\mathrm{J}}}^{\mathrm{d}}>0$
are the disease costs for adults and juveniles, respectively, and
$c^{\mathrm{v}}>0$ is the vaccination cost. We define the relative costs
$\varepsilon=c_{{\mathrm{J}}}^{\mathrm{d}}/c_{{\mathrm{A}}}^{\mathrm{d}}$ and
$r=c^{\mathrm{v}}/c_{{\mathrm{A}}}^{\mathrm{d}}$. Upon normalization, we will
assume from now on that $c_{{\mathrm{A}}}^{\mathrm{d}}=1$. The fraction of the
vaccination cost supported by the society is given by $\delta\in[0,1]$, where
$\delta=1$ means that all cost is supported by the society (normally, the
State), where $\delta=0$ means that the entire cost of the vaccination is paid
by the vaccinated individual.

Note that $I^{\mathrm{en}}_{{{\mathrm{A}}},{{\mathrm{J}}}}$ depend explicitly
on $p$ and $\lambda$, cf. Thm. 1.

Parameter | Description
---|---
$c_{{\mathrm{A}}}^{\mathrm{d}}$ | disease cost of an adult
$c_{{\mathrm{J}}}^{\mathrm{d}}$ | disease cost of a juvenile
$c^{\mathrm{v}}$ | vaccination cost
$\delta$ | Fraction of the vaccination costs supported by the society
$\varepsilon:=c_{{\mathrm{J}}}^{\mathrm{d}}/c_{{\mathrm{A}}}^{\mathrm{d}}$ | relative disease cost of juveniles vs. adults
$r:=c^{\mathrm{v}}/c_{{\mathrm{A}}}^{\mathrm{d}}$ | relative vaccination cost vs. adults disease cost
Table 3: Cost variables used in the model. Upon normalization
$c_{{\mathrm{A}}}^{\mathrm{d}}=1$, results presented in this article will
depend only on $\delta$, a modeling parameter, $\varepsilon$ and $r$. The
values for the relative costs $\varepsilon$ and $r$ used in this work are
arbitrary and used for illustration purposes.

We define the acceptable social-cost region as

$\mathcal{V}=\left\\{(p,\lambda)\in[0,1]\times[0,1]:\Phi_{\varepsilon,r,\delta}(p,\lambda):=\phi(p,\lambda)-\phi_{0}\leq
0\right\\}\ ,$

where $\phi_{0}=\phi(0,0)$ is the social cost of the disease in an
unvaccinated population.

We define two critical values: $\lambda_{\sup}$, below which social-cost
acceptance depends on vaccine coverage $p$; and $\lambda_{\inf}$, below which
social-cost is unacceptable for any vaccine coverage $p$.

$\lambda_{\sup}=\sup_{\Phi(p,\lambda)>0}\lambda,\quad\lambda_{\inf}=\inf_{\Phi(p,\lambda)<0}\lambda\
.$

Figure 2: The light-blue region indicates the acceptable cost region $\Phi<0$,
while the grey region is the disease-free region ${\mathcal{R}}_{p}<1$. The
number $\lambda_{\inf}$ indicates the minimum value of vaccine efficacy such
that a sufficiently high vaccine coverage will guarantee that the disease has
an acceptable social cost at equilibrium. The number $\lambda_{\sup}$
indicates the minimum value of $\lambda$ such that _any_ vaccine coverage is
in the acceptable social-cost region. We assume a juvenile/adult relative cost
$\varepsilon=0.15$, a vaccine/disease cost $r=0.1$, and all vaccination costs
are supported by the vaccinated individual, i.e., $\delta=0$.

Fig. 2 illustrates the acceptable social-cost region in the parameter space
$(p,\lambda)$ when $\delta=0$. Note that there is a subregion in which is is
possible to eliminate the disease, i.e., ${\mathcal{R}}_{p}<1$.

###  2.3 Individual cost and Nash equilibria

Following [2], we assume that individuals freely choose to be vaccinated
according to the perceived relative costs of the disease and of the
vaccination. For each $(p,\lambda)$, let us define
$\Pi^{{\mathrm{nv}}}_{{{\mathrm{A}}}}$ and
$\Pi^{{\mathrm{nv}}}_{{{\mathrm{J}}}}$ as the stationary (i.e., at
equilibrium) probabilities of getting the disease as an adult and as a
juvenile for unvaccinated individuals; and
$\Pi^{{\mathrm{v}}}_{{{\mathrm{A}}}}$ to be the stationary probability of
getting the disease as an adult, if vaccinated at birth. These values are
equal to zero at the disease-free equilibrium and non-zero at the endemic
equilibrium. Furthermore, they are continuous functions from the model
parameters, cf. [18, Sec. 3.4].

We obtain explicit expressions for each of these three parameters.

For $\Pi_{{\mathrm{J}}}^{{\mathrm{nv}}}$, we consider a given individual in
the class $S_{{\mathrm{J}}}$, from which there are two possible exits. Either
that given individual contracts the disease (and move to the class
$I_{{\mathrm{J}}}$) or he or she turns into an adult without being infected
and moves to the class $S_{{\mathrm{A}}}$. Explicitly,

$\Pi_{{\mathrm{J}}}^{{\mathrm{nv}}}(p,\lambda):=\frac{\beta
I^{{{}^{*}}}S^{{{}^{*}}}_{{\mathrm{J}}}}{(\beta
I^{{{}^{*}}}+\nu)S^{{{}^{*}}}_{{\mathrm{J}}}}=\frac{\beta I^{{{}^{*}}}}{\beta
I^{{{}^{*}}}+\nu}\ .$

The probability that a non-vaccinated adult gets the disease is given by the
probability that a previously non-vaccinated juvenile does not get the disease
as a juvenile and then gets the disease as an adult. Therefore

$\Pi_{{\mathrm{A}}}^{{\mathrm{nv}}}:=\left(1-\Pi_{{\mathrm{J}}}^{{\mathrm{nv}}}\right)\frac{\beta
I^{{{}^{*}}}S^{{{}^{*}}}_{{\mathrm{A}}}}{(\beta
I^{{{}^{*}}}+\mu)S^{{{}^{*}}}_{{\mathrm{A}}}}=\frac{\nu}{\beta
I^{{{}^{*}}}+\nu}\frac{\beta I^{{{}^{*}}}}{\beta I^{{{}^{*}}}+\mu}\ ,$

with $I^{*}:=I^{*}_{{\mathrm{J}}}+I^{*}_{{\mathrm{A}}}$. Finally, the
probability that a vaccinated adult gets the disease is the probability that
the vaccine is effective only during the juvenile phase, $1-\lambda$, times
the probability to get the disease from the class $S_{{\mathrm{A}}}$, i.e.,

$\Pi_{{\mathrm{A}}}^{{\mathrm{v}}}:=(1-\lambda)\frac{\beta I^{{{}^{*}}}}{\beta
I^{{{}^{*}}}+\mu}\ .$

We define the _individual cost function_ at endemic equilibrium, which
corresponds to the expected cost of the individual strategy of being
vaccinating with probability $q$ in a population with coverage $p$:

$\displaystyle\Psi_{\varepsilon,r,\delta}(q,p,\lambda)$
$\displaystyle:=(1-q)(\Pi^{{\mathrm{nv}}}_{{{\mathrm{A}}}}+\varepsilon\Pi^{{\mathrm{nv}}}_{{\mathrm{J}}})+q(\Pi^{\mathrm{v}}_{{\mathrm{A}}}+r(1-\delta))$
$\displaystyle=\Pi^{{\mathrm{nv}}}_{{{\mathrm{A}}}}+\varepsilon\Pi^{{\mathrm{nv}}}_{{\mathrm{J}}}+q\left[r(1-\delta)-\pi(p,\lambda)\right]\
,$

where the _vaccination-infection risk index_ , introduced in [12], is given by

$\pi(p,\lambda):=\Pi^{{\mathrm{nv}}}_{{{\mathrm{A}}}}(p,\lambda)+\varepsilon\Pi^{{\mathrm{nv}}}_{{\mathrm{J}}}(p,\lambda)-\Pi^{\mathrm{v}}_{{\mathrm{A}}}(p,\lambda)\
.$

The individual _vaccination marginal expected payoff gain_ $E(q,p)$ of an
individual that uses the strategy of vaccinating with probability $q$ in a
population that vaccinates with probability $p$ is given by

$E(q,p):=E(q,p;\varepsilon,r,\delta,\lambda):=\Psi_{\varepsilon,r,\delta}(0,p,\lambda)-\Psi_{\varepsilon,r,\delta}(q,p,\lambda)\
.$

######  Definition 1.

The population vaccination strategy $p_{*}$ is a vaccination _Nash
equilibrium_ , if

$\displaystyle
E(q,p_{*})-E(p_{*},p_{*})=(p_{*}-q)\left[r(1-\delta)-\pi(p_{*},\lambda)\right]\leq
0,$

for every strategy $q\in[0,1]$.

In simple words, we say that the system is at Nash equilibrium if the
vaccination coverage $p_{*}$ is such that for every individual that uses a
strategy $q$ the expected payoff is not larger than the one it would have if
the strategy $p_{*}$ were used.

######  Proposition 4.

The model given by Eqs. (2)–(7) has at least one Nash equilibrium.

###### Proof.

If $\pi(0,\lambda)\leq r(1-\delta)$, then $p_{*}=0$ is a Nash equilibrium. If
$\pi(1,\lambda)\geq r(1-\delta)$, then $p_{*}=1$ is a Nash equilibrium. If
both inequalities are false there is at least one value of $p_{*}\in(0,1)$
such that $\pi(p_{*},\lambda)=r(1-\delta)$ and $p_{*}$ is a Nash equilibrium.
∎

For high vaccine efficacy $\lambda>\lambda^{*}$ and $\delta\in[0,1)$, the
vaccination coverage that results from individuals’ choices is below the
elimination threshold $p_{c}$, defined in Prop. 3.

######  Proposition 5.

Let $\varepsilon,r>0$, $\delta\in[0,1)$, $\lambda\in[\lambda_{\mathrm{c}},1]$,
where $\lambda_{\mathrm{c}}$ is given by Prop. 3. Let
$p_{\mathrm{c}}^{\lambda}$ given by Prop. 3 and $p_{*}^{\lambda}$ a Nash
equilibrium of the associated model. Then,
$p_{*}^{\lambda}<p_{\mathrm{c}}^{\lambda}$.

###### Proof.

From Prop. 3, for any value $p>p_{\mathrm{c}}^{\lambda}$ it is true that
$\pi(p,\lambda)=0$. From the continuity of $\pi$, we conclude that
$\pi(p_{c}^{\lambda},\lambda)=0$. Assume that $p_{*}^{\lambda}\geq
p_{\mathrm{c}}^{\lambda}>0$. From Def. 1 we have that
$(p_{*}^{\lambda}-q)\left[r(1-\delta)-\pi(p_{*}^{\lambda},\lambda)\right]\leq
0$ for every $q\in[0,1]$, therefore $p_{*}^{\lambda}\leq q$, for every
$q\in[0,1]$, which is a contradiction. ∎

Note that this result generalizes for arbitrary efficacy $\lambda$ the idea,
already present in [2], that a Nash equilibrium of a vaccination game is
always below the threshold to eradicate a disease.

Inspired by the concept of evolutionary stable strategy in game dynamics, cf.
[7], we define:

######  Definition 2.

The population vaccination strategy $p_{*}$ is an _evolutionary stable
vaccination (ESV)_ strategy, if there is a $\tau_{0}>0$, such that for every
$\tau\in(0,\tau_{0})$ and for every $q\in[0,1]$, with $q\neq p_{*}$,

$E(q,(1-\tau)p_{*}+\tau q)-E(p_{*},(1-\tau)p_{*}+\tau q)<0.$

We are ready to state the conditions for the Nash equilibrium to be ESV.

######  Proposition 6.

Let $p_{*}$ be a Nash equilibrium of the vaccination game. If $p_{*}=0$ or
$p_{*}=1$, then $p_{*}$ is an ESV. Furthermore, if
$\pi(p_{*},\lambda)=r(1-\delta)$, $p_{*}$ is an ESV if and only if
$\pi(p,\lambda)$ is decreasing at $p=p_{*}$. In particular, $p_{*}\in(0,1)$ is
an ESV if and only if $\pi(p,\lambda)$ is decreasing at $p=p_{*}$.

###### Proof.

This proof follows ideas from [12]. Let $p_{*}=0$ ($p_{*}=1$) be a Nash
equilibrium. From Def. (1), we conclude that $\pi(0,\lambda)\leq r(1-\delta)$
($\pi(1,\lambda)\geq r(1-\delta)$, respect.). Assume that a strict inequality
is valid. Let $\tau_{0}$ be small enough such that for all $\tau<\tau_{0}$, it
is valid that $\pi(\tau q,\lambda)<r(1-\delta)$
($\pi(1-\tau(1-q),\lambda)>r(1-\delta)$, respect.). It is clear that $E(q,\tau
q)-E(0,\tau q)=-q(r(1-\delta)-\pi(\tau q,\lambda))<0$
($E(q,1-\tau(1-q))-E(1,1-\tau(1-q))=(1-q)(r(1-\delta)-\pi(1-\tau(1-q),\lambda))<0$,
respect.), for all $q\neq p_{*}$.

For the second part, note that $\pi(p_{*},\lambda)=r(1-\delta)$ implies that

$\displaystyle E(q,(1-\tau)p_{*}+\tau q)-$ $\displaystyle
E(p_{*},(1-\tau)p_{*}+\tau q)$
$\displaystyle=-(q-p_{*})(\pi(p_{*},\lambda)-\pi((1-\tau)p_{*}+\tau
q,\lambda))\ ,$

and therefore $p_{*}$ is an ESV if and only if $\pi$ is decreasing in the
first argument at $p=p_{*}$. Finally, the final result follows from Def. 1. ∎

(a)

(b)

Figure 3: Nash equilibria as a function of vaccine efficacy $\lambda$ and
vaccination coverage $p$ for relative vaccination costs (a) $r=0.25$ and (b)
$r=0.30$. The Light-red region is such that $\pi(p)>r(1-\delta)$, i.e., in
this region a rational individual will accept to be vaccinated with a
probability larger than the population average. In particular, in that region,
there is an individual incentive to increase the overall vaccination coverage.
Red dashed and full lines correspond to unstable and stable Nash equilibria,
respectively. The function $\pi$ is decreasing (increasing) with respect to
$p$ in the full (dashed) red line, cf. Prop. 6. The grey region is the
disease-free region where ${\mathcal{R}}_{p}<1$, and the full black line is
the disease-free threshold ${\mathcal{R}}_{p}=1$. The horizontal black dotted
line exemplifies the dynamics of rational individuals (indicated by the
arrows) assuming a vaccine efficacy of $\lambda=0.70$. The region between
$\lambda_{\inf}^{\mathrm{bi}}$ and $\lambda_{\sup}^{\mathrm{bi}}$ is the
region for model bistability, in which we find three Nash equilibria, two
stable and one unstable in between. We assume a juvenile/adult relative cost
$\varepsilon=0.15$ and all vaccination costs are supported by the vaccinated
individual $\delta=0$. Note that it is not possible to reach the disease-free
region through voluntary vaccination if there is no incentive to be
vaccinated. However, in case (a), the region in which there is no individual
incentive to increase the vaccination coverage close to the disease-free
region is disconnected from the set of vaccination coverage $p=0$.

Fig. 3 illustrates two possible situations described in Prop. 6: (a) the two
pure strategies are ESV and there exists an interior Nash equilibrium that is
unstable; (b) for higher relative vaccination costs the interior Nash
equilibrium is stable when condition on Prop. 6 is met. For both situations
described, there is a range
$\left(\lambda_{\inf}^{\mathrm{bi}},\lambda_{\sup}^{\mathrm{bi}}\right)$ for
the vaccine efficacy $\lambda$ were the model presents bi-stability.

##  3 Discussion and Numerical Examples

In this section, we present several numerical examples to discuss the present
work. Parameters will be, except otherwise said, taken from Table 1. In
particular, chickenpox epidemiology fits our framework, as it is a mild
disease for children that can have increased risk for adults and its use in a
universal vaccination program is debatable [9]. However, the framework
developed here may be applied to several other situations, such as Zika or
rubella.

Figure 4: Proportion of infectious individuals at endemic equilibrium without
vaccination ($p=0$) as a function of the basic reproduction number
${\mathcal{R}}_{0}$. The full line indicates
$I^{\mathrm{en}}=I^{\mathrm{en}}_{{\mathrm{J}}}+I^{\mathrm{en}}_{{\mathrm{A}}}$,
while the dotted and dashed lines indicate $I^{\mathrm{en}}_{{\mathrm{J}}}$
and $I^{\mathrm{en}}_{{\mathrm{A}}}$. Note that for larger values of
${\mathcal{R}}_{0}$, the fraction of juveniles approaches the full number of
infectious, indicating that a highly-transmissible disease with permanent
immunity will be, in the stationary state, a childhood disease. However, when
the transmission is low, a significant number of infectious individuals is
adult.

Fig. 4 shows the proportion of infectious individuals at equilibrium as a
function of the basic reproduction number without vaccination, i.e., $p=0$.
The total proportion of infectious individuals $I^{\mathrm{en}}$ in the
endemic equilibrium is an increasing function, as is the case of the
proportion of juveniles $I_{{\mathrm{J}}}^{\mathrm{en}}$ in the same
equilibrium. The fraction of adults increases for small values of
${\mathcal{R}}_{0}$ and then decreases. We conclude that a highly
transmissible disease associated with permanent immunity will be, in
equilibrium, a childhood disease. If the effect of this disease is mild in
juveniles, there is no severe economic cost associated with the endemic state.
This is the main reason we will always compare the economic cost associated
with a vaccine program with vaccine efficacy $\lambda$ and vaccine coverage
$p$ with the no-vaccination endemic state, cf. definition of $\Phi$ in
Subsection 2.2. For that choice of parameters, most of the infectious
individuals are below 15 years old, but a reasonable proportion of infectious
individuals is above this value.

(a)

(b)

Figure 5: (a) Fraction of infectious individuals at endemic equilibrium
assuming vaccination coverage $p$ and vaccine efficacy $\lambda$. (b) Fraction
of infected juveniles, with respect to the number of infected individuals at
endemic equilibrium, $I_{{\mathrm{J}}}^{{\mathrm{en}}}/I^{{\mathrm{en}}}$, as
a function of the vaccine coverage $p$ and vaccine efficacy $\lambda$. Note
that increasing the vaccine coverage implies a smaller number of infected
individuals but the disease became more relevant among adults. The grey region
in the upper left corner of both graphs indicates the disease-free region.

Assuming ${\mathcal{R}}_{0}=8$, the inclusion of a vaccination scheme is
illustrated in Fig. 5. It clearly shows that for the relevant set of
parameters, the inclusion of a vaccination program will decrease the overall
number of infectious individuals in the endemic equilibrium but it will
increase the fraction of adults. Therefore, the introduction of the
vaccination scheme should be pondered to avoid negative outcomes for the
population.

After the introduction of the vaccination, two natural questions arise: i)
_are people willing to be vaccinated?_ , and ii) _has the individual behavior
a positive or negative effect on society?_ The first question is addressed by
introducing an individual cost of being vaccinated (that includes the
perceived risk of the vaccine, eventual absence to work to go or to take the
children to the vaccination site, the financial cost of buying the vaccine,
etc) and the cost of non-being vaccinated, i.e., all the costs associated to
contracting the disease. If the first is larger, then rational individuals
will be vaccinated, if it is smaller, they will not. The equality points
correspond to the Nash equilibrium of the model. For the second question, we
discuss if a given vaccination is in the acceptable social cost region.
Ideally, we shall try to find a stable Nash equilibrium within that region,
i.e., with $\Phi<0$. However, this is not always possible.

(a)

(b)

(c)

Figure 6:  Social _vs._ individual interest with vaccine coverage $p$ and
vaccine efficacy $\lambda$, with fixed $\varepsilon=0.15$. The grey region
marks the disease-free region, i.e. ${\mathcal{R}}_{p}<1$. We consider
different scenarios: (a) with all vaccination costs assumed by the individual,
i.e., $\delta=0$, with a high-cost vaccine $r=0.35$; (b) $\delta=0$ and low-
cost vaccine $r=0.01$; or (c) with shared costs between the individual and the
society $\delta=0.36$ for the high-cost vaccine $r=0.35$. In light-blue
region, light-purple region and grey region the social cost is lower than the
social cost of an unvaccinated population, i.e., $\Phi<0$; the blue continuous
line indicates the level set $\Phi=0$. In the light-red region and light-
purple region, it is in the individual interest to be vaccinated with a larger
probability than the population average; in the light-red region the social
cost of vaccination is positive. The red line corresponds to the set of Nash
equilibria. The horizontal red dotted lines represent the range
$\lambda\in\left(\lambda_{\inf}^{\mathrm{bi}},\lambda_{\sup}^{\mathrm{bi}}\right)$
where the model has bi-stability and the horizontal blue dotted lines
represent the range $\lambda\in\left(\lambda_{\inf},\lambda_{\sup}\right)$
were the social cost of vaccination is acceptable only if the level of
vaccination is sufficiently high. In (c) the vaccine costs are shared so that
$\lambda_{\sup}^{\mathrm{bi}}=\lambda_{\sup}$ and the blue and red lines
intersect for $p=0$. In particular, all Nash equilibria $p>0$ are in the
socially cost-efficient region.

Figs. 6(a) and 6(b) illustrate the regions on the parameter space
$(p,\lambda)$ in which the individual and the social interests coincide and
differ when all the individual vaccination costs are assumed by the
beneficiary (i.e., there are no shared costs, as, for example, government
subsidies). Fig. 6(c) introduces shared costs for high-cost vaccines, i.e.,
the society absorbs part of the individual cost.

Close to the disease-free region ${\mathcal{R}}_{p}<1$, there is always a
barrier where there is no individual interest to be vaccinated, as the
infection rates at that region will be residual. In Fig. 6(a), the region
where there is a social interest in increasing the vaccination, but there is
an individual rejection of it, is a connected set. In Fig. 6(c) this region is
disconnected. In the former case, it is possible that a decrease in the value
of $\lambda$ (something not included in our model, but that may happen due,
for instance, to the introduction of new variants or simply because it’s
perceived as so by the population) causes a near-perfect vaccination scheme to
collapse into a non-vaccination situation (i.e., with $p\approx 0$), due only
to rational individual behavior. This is not possible when this region becomes
disconnected, bringing extra stability to a near-optimal vaccination scheme.

In Fig. 6 we indicate the values of $\lambda_{\inf}$ ($\lambda_{\sup}$), the
minimum efficacy such that there is an individual incentive to be vaccinated
for $p$ large enough (for any value of $p$, respect.) and
$\lambda_{\inf,\sup}^{\mathrm{bi}}$ the minimum and maximum values to the
existence of bi-stable Nash vaccination equilibrium. Depending on the costs of
the vaccine for society and for individuals, their interests may not always
agree: for a certain range of vaccine efficacy, it may be favorable for
society to increase vaccination coverage, but due to the high cost of the
vaccine, individuals choose not to be vaccinated, cf. Fig. 6(a) for
$\lambda\in(\lambda_{\inf},\lambda_{\inf}^{\mathrm{bi}})$. For a different set
of parameters it may be favorable for individuals to vaccinate, due to the low
vaccine cost, but not be beneficial for society, cf. Fig. 6(b) for
$\lambda\in(\lambda_{\inf}^{\mathrm{bi}},\lambda_{\inf})$. This situation can
be changed by allowing the vaccination costs to be shared. For example, in
Fig. 6(c), $\delta$ was chosen such that
$\lambda_{\sup}^{\mathrm{bi}}=\lambda_{\sup}$. In this case, individual
vaccination is enhanced for lower vaccine efficacy, as compared to Fig. 6(a).
Moreover, all Nash equilibrium $p>0$ are in the acceptable social cost region,
avoiding individuals choosing to be vaccinated where their choice would
increase social costs.

Figure 7: Consider the space of parameter given by $\delta$, the fraction of
the vaccination cost supported by the society, and $p$, the vaccination
coverage. The light-blue and light-purple are the regions such that the level
of vaccination has a social cost lower than of no vaccination, while the
light-red and light-purple regions are such that a rational individual will
choose to be vaccinated with a larger probability than the average individual.
The light-purple region is the objective of the health authorities, where
individuals freely decide to be vaccinated and the overall coverage is cost-
efficient, i.e., $\Phi<0$. We assume $\varepsilon=0.15$ and $\lambda=0.6$.
Consider the example illustrated by the yellow arrows. If we start in (a) with
$p=0.5$ (indicated by a dotted vertical line) and $\delta=0$, rational
individuals will not vaccinate, but vaccination will benefit society. However,
if the vaccination cost starts to be shared between society and beneficiary,
increasing the value of $\delta$ to above approx 0.5 [point (b)], rational
individuals will start to vaccinate. Vaccination coverage will increase until
close to 1, inside the acceptable socially-cost region and the shared costs
can be relaxed to $\delta=0.3$ [point (c)]. The chosen yellow points
correspond to Fig. 8(a), 8(b) and 8(c), respectively.

(a)

(b)

(c)

Figure 8: Parameter space $(p,\lambda)$ for different values of the fraction
of the cost supported by the society: (a) $\delta=0$, (b) $\delta=0.5$, (c)
$\delta=0.3$. In light-blue and light-purple indicates the region where the
level of vaccination has a social cost lower than of no vaccination, while
light-red and light-purple indicates the region in which rational individuals
find incentives to have themselves vaccinated. Neutral social-cost threshold
$\Phi=0$ and Nash equilibria are indicated by blue and red curves,
respectively. Red solid lines indicate stable Nash equilibria, while dotted
lines indicated unstable equilibrium. The vertical black dotted line
represents the level of vaccination $p=0.5$ in (a) and (b) and $p=1$ in (c).
The horizontal black dotted line represents the vaccine efficacy
$\lambda=0.6$, corresponding to the examples studied in points (a), (b), and
(c) in Fig. 7.

The effects of sharing costs are further explored in Figs. 7 and 8. Fig. 7
shows a particular example, highlighted by the yellow arrows. Starting with a
vaccination coverage of $p=0.5$, by changing the value of $\delta$, it is
possible to create incentives such that rational individuals accept to be
vaccinated, moving from the light-blue region to the light-red region, i.e.
from point (a) to point (b). The natural dynamics will lead the population to
a state in which the level of vaccination is high and the population is in a
cost-efficient equilibrium. After that, it is possible to decrease vaccination
incentives without decreasing vaccination coverage, i.e., moving towards point
(c). Fig. 8(a), 8(b) and 8(c), show the superposition of the individual and
social interests for different scenarios for shared costs corresponding to the
three points depicted in Fig. 7, respectively.

##  4 Summary of conclusions

This work starts from the fact that imperfect vaccination can be worst than no
vaccination for a specific group of diseases to discuss the implementation of
specific strategies that induce rational individuals to be vaccinated in a
socially cost-efficient way. This is an important issue, as, at least in
developed democracies, forced vaccination is considered unacceptable, but
positive and negative incentives to boost vaccination coverage are routinely
used. In fact, many restrictions to non-vaccinated individuals were
implemented during the COVID-19 pandemic, even in developed democracies,
showing that these strategies are always been considered, at least in extreme
cases.

For many reasons, in particular, the almost-extinction of many vaccine-
preventable diseases, vaccine-skeptical groups are present in almost every
country. In this work, we introduced an age-structured model, considered
different effects of the disease in adults and juveniles, considered imperfect
vaccines, and study socially efficient vaccine coverage, Nash equilibrium
vaccination strategies, and, more importantly, the intersection between these
two groups. Finally, we show how sharing the cost between individuals and
society can boost vaccination coverage, moving, if not to the extinction of
the disease, at least to its long-term control.

## Acknowledgements

This work is funded by national funds through the FCT – Fundação para a
Ciência e a Tecnologia, I.P., under the scope of the projects UIDB/00297/2020
and UIDP/00297/2020 (Center for Mathematics and Applications – NOVA Math)
and 2022.03091.PTDC _Mathematical Modelling of Multi-scale Control Systems:
applications to human diseases_ (CoSysM3).

## References

  * [1] A. I. Abushouk, A. Negida, and H. Ahmed.  An updated review of zika virus.  J. Clin. Virol., 84:53–58, NOV 2016.
  * [2] C. Bauch and D. Earn.  Vaccination and the theory of games.  PNAS, 101(36):13391–13394, 2004.
  * [3] S. L. Chang, M. Piraveenan, P. Pattison, and M. Prokopenko.  Game theoretic modelling of infectious disease dynamics and intervention methods: a review.  J. Biol. Dynam., 14(1):57–89, jan 2020.
  * [4] A. J. Daley, S. Thorpe, and S. M. Garland.  Varicella and the pregnant woman: Prevention and management.  Aust. NZ J. Obstetrics & Gynaecology, 48(1):26–33, FEB 2008.
  * [5] P. Doutor, P. Rodrigues, M. C. Soares, and F. A. C. C. Chalub.  Optimal vaccination strategies and rational behaviour in seasonal epidemics.  J. Math. Biol., 73:1437–1465, 2016.
  * [6] P. FINE.  Herd-immunity - history, theory, practice.  Epidemiol. Rev., 15(2):265–302, 1993.
  * [7] J. Hofbauer and K. Sigmund.  Evolutionary Games and Population Dynamics.  Cambridge University Press, 1998.
  * [8] L. Laguzet.  High order variational numerical schemes with application to Nash-MFG vaccination games.  Ric. Mat., 67(1):247–269, 2018.
  * [9] Y. H. Lee, Y. J. Choe, J. Lee, E. Kim, J. Y. Lee, K. Hong, Y. Yoon, and Y.-K. Kim.  Global varicella vaccination programs.  Clin. Exp. Pediatr., 65(12):555–562, DEC 2022.
  * [10] J. Liu, B. F. Kochin, Y. I. Tekle, and A. P. Galvani.  Epidemiological game-theory dynamics of chickenpox vaccination in the USA and Israel.  J. R. Soc. Interface, 9(66):68–76, jan 2012.
  * [11] P. Manfredi and A. D’Onofrio.  Modeling the interplay between human behavior and the spread of infectious diseases.  Springer New York, Oct. 2013.
  * [12] J. Martins and A. Pinto.  Bistability of Evolutionary Stable Vaccination Strategies in the Reinfection SIRI Model.  B. Math. Biol., 79:853–883, 2017.
  * [13] T. Panagiotopoulos, I. Antoniadou, E. Valassi-Adam, and A. Berger.  Increase in congenital rubella occurrence after immunisation in Greece: retrospective survey and systematic review How does herd immunity work?  BMJ, 319(7223):1462–1467, 1999.
  * [14] M. Riera-Montes, K. Bollaerts, U. Heininger, N. Hens, G. Gabutti, A. Gil, B. Nozad, G. Mirinaviciute, E. Flem, A. Souverain, T. Verstraeten, and S. Hartwig.  Estimation of the burden of varicella in europe before the introduction of universal childhood immunization.  BMC Infect. Dis., 17, MAY 18 2017.
  * [15] F. Salvarani and G. Turinici.  Optimal individual strategies for influenza vaccines with imperfect efficacy and durability of protection.  Math. Biosci. Eng., 15(3):629–652, 2018.
  * [16] Z. Shuai and P. van den Driessche.  Global stability of infectious disease models using lyapunov functions.  SIAM J. Appl. Math., 73(4):1513–1532, 2013.
  * [17] P. van den Driessche and J. Watmough.  Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission.  Math. Biosci., 180(SI):29–48, 2002.
  * [18] M. Viana and J. M. Espinar.  Differential equations—a dynamical systems approach to theory and practice, volume 212 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, 2021.
  * [19] Z. Wang, C. T. Bauch, S. Bhattacharyya, A. d’Onofrio, P. Manfredi, M. Perc, N. Perra, M. Salathé, and D. Zhao.  Statistical physics of vaccination.  Physics Reports, 664:1 – 113, 2016.