# Exploring Music Genre Classification: Algorithm Analysis and Deployment
Architecture

1st Ayan Biswas1, 2nd Supriya Dhabal2, and 3rd Palaniandavar Venkateswaran3
2Electronics & Communication Engineering Department,
Netaji Subhash Engineering College, Kolkata - 700152, INDIA
<EMAIL_ADDRESS>13Dept. of Electronics and Tele-communication
Engineering ,
Jadavpur University, Kolkata – 700 032, INDIA
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Music genre classification has become increasingly critical with the advent of
various streaming applications. Nowadays, we find it impossible to imagine
using the artist’s name and song title to search for music in a sophisticated
music app. It is always difficult to classify music correctly because the
information linked to music, such as region, artist, album, or non-album, is
so variable. This paper presents a study on music genre classification using a
combination of Digital Signal Processing (DSP) and Deep Learning (DL)
techniques. A novel algorithm is proposed that utilizes both DSP and DL
methods to extract relevant features from audio signals and classify them into
various genres. The algorithm was tested on the GTZAN dataset and achieved
high accuracy. An end-to-end deployment architecture is also proposed for
integration into music-related applications. The performance of the algorithm
is analyzed and future directions for improvement are discussed. The proposed
DSP and DL-based music genre classification algorithm and deployment
architecture demonstrate a promising approach for music genre classification.

###### Index Terms:

Music Genre Classification, Feature Extraction, Deep Learning

##  I Introduction

Aside from providing entertainment, music is one of the easiest ways to
communicate among people, a way to share emotions, and a place to keep
memories and emotions. Emotions can be expressed succinctly and effectively
through music. Depending on the mood and objective of the listener, people
select different music at different times. As Internet technology flourishes,
more and more music is available on personal computers, in music libraries,
and via the Internet. Systems that can automatically analyze music, like
categorizing it, searching through it, and creating playlists, are crucial for
efficiently managing music. Several studies have proposed that music mood can
also be used to classify and recommend music. There are a number of existing
mood-based music recommendation systems that categorize some moods and map
those moods into discrete regions in two or three dimensions. We have also
included a section that shows a possible architecture for deploying the
solution to mobile applications or the web since as developers of music
applications it is pretty ambiguous whether these scientific solutions should
be deployed or not.

The remainder of this paper is organized as follows: Section II reviews the
related works. Section III presents the overall approach, the deep-learning
model, the algorithm and the discussion on the final results. The deployment
architecture of the music genre classification algorithm has been discussed in
Section IV. Finally, Section V draws the concluding remarks of our paper.

##  II Related Works

Recently, there has been an increase in the attention given to analyzing audio
to extract different kinds of information, specifically in relation to music
and emotions. Research has focused on developing automated methods for
classifying music according to its mood or emotional content. Some different
proposed approaches are there, including the use of spectral and harmonic
features to infer the mood of a music piece. These features have been linked
to human perception of music and moods, and have been used to classify music
according to different mood labels using neural networks. This literature
review suggests that the use of spectral and harmonic features, along with
neural network-based classification methods, can be a promising approach for
classifying music according to mood.

Bhat et al. have proposed a number of different approaches to solving the
problem in their work [1], including using spectral and harmonic features to
infer the mood of a given music piece. In particular, features such as rhythm,
harmony, spectral feature, and others have been studied in order to classify
songs according to their mood. This has been based on Thayer’s model, which
proposes that certain features of music are linked to human perception of
music and moods.

In this paper [2] by Kim et al., a probability-based music mood model and its
application to a music recommendation system have been presented. In this
approach three types of mood-based music recommendation players, for PC and
mobile devices, and the web has been implemented. This paper also shows the
analysis result of users’ satisfaction and mood reappearance test after
listening to music.

According to this paper [3] by Patel et al., sound is the most important
aspect of this project and can be distinguished by its pitch, quality, and
loudness. The fundamental tone and the harmonics are generated and give rise
to different musical notes. The Fourier transform has been used for breaking
musical tones into sinusoidal waves.

Tzanetakis et al. in [4] demonstrated that music genre classification can be
done by manipulating three types of features that represent the texture,
rhythm, and pitch of the music. They evaluated the effectiveness and
significance of these features by training machine learning models using real-
world audio collections in their research.

In prior literature, [3], the development of an algorithm for the
identification of musical notes was presented, yet the crucial aspect of
deployment architecture remained elusive. This study aims to bridge that gap
by proffering a comprehensive deployment schema that ensures optimal
performance and minimal error in the identification of musical notes. The
succeeding sections of this paper are dedicated to delving into the
intricacies of the proposed implementation, providing a detailed account of
its deployment and execution.

##  III Our Approach

###  III-A Methodology

The music signal feature classification begins by recording a music sound and
obtaining the corresponding waveform [5]. The frequency of the notes within
the music is identified by analyzing the duration of each note in the time
domain. An averaging process is applied to reduce the number of samples and
fluctuations. The envelope of the original signal can also be extracted.
Subsequently, thresholding is performed to establish a threshold value for
identifying the maximum peaks in the signal. A technique of dynamic adaptive
thresholding [6], which adjusts the threshold value based on the number of
peaks, can be used for this purpose. Next, a width interval is selected to
facilitate further operations. The width interval is chosen such that a larger
number of peaks can be condensed within a smaller length. The width interval
increases as the sampling frequency increases, as it is an essential aspect of
the sampling process. To find the sine waves in a signal, a Fourier transform
is applied, and zero padding is used to minimize error and get the Discrete
Fourier Transform (DFT) of the signal. The frequency of the musical notes is
identified by analyzing the frequency of the resulting signal from the DFT.

###  III-B Feature Extraction from Music Samples

The focus on extracting as many features [7, 8, 9] as possible is motivated by
the fact that this can make the subsequent classification task more
straightforward. Some common features that are extracted from audio signals
include tonality, pitch, temporal energy, harmonicity, spectral centroid, and
Mel-Frequency Cepstral Coefficients (MFCC) [10].

Figure 1: Architecture of the Feature Extractor and Classifier

There are two main types of features that describe an audio signal: global
descriptors and instantaneous descriptors. Global descriptors are computed for
the entire signal and help identify steady patterns in the signal, such as the
total energy of an audio clip or the emotional tone of a song.

Instantaneous descriptors provide information about the dynamic and temporal
variations of a signal. These descriptors are obtained by dividing the signal
into small segments, called frames, and then applying pre-processing
techniques to each frame. The features calculated for each frame are usually
related to time, spectral shape, harmonic, and energy. This paper focuses on
extracting instantaneous descriptors for each frame, discussed in the work of
[11] and [10].

####  III-B1 Pitch

Pitch is a metric that describes the regularity of a sound wave or the
perceived fundamental frequency of the signal. The true frequency of the
signal can be determined precisely, but it may not match the perceived pitch
due to the presence of harmonics. To determine the pitch, the auto-correlation
sequence (ACS) for a given frame of the signal is calculated using a specific
formula as per equation (1).

Figure 2: Plot of the pitch for all genres of music samples

$r(m)=\frac{1}{N}\sum_{n=0}^{N-|m|-1}x(n+|m|)x(n)$ (1)

where $N$ is the length of the frame in samples and $x$ is the input signal,
such as speech or audio signal

####  III-B2 Temporal Energy

The temporal energy E, which is a measure of the strength of the signal over a
specific frame of time, is calculated by finding the average of the squared
values of the signal over that frame. This is expressed mathematically in
equation (2). The energy feature can be used to distinguish between voiced
frames, which contain significant information about the signal, and unvoiced
frames, which are typically silent or noise-like, by comparing the energy
values to a fixed threshold value.

$E=\frac{1}{N}\sum_{n=0}^{N-1}x^{2}(n)$ (2)

Figure 3: Temporal Energy measurement of a random music signal from GTZAN
dataset

####  III-B3 Tonality Measure

A significant amount of background noise or sensor noise can obscure the true
tone of an audio or speech signal. Tonality is a metric that describes how
much of the signal has a tone-like or noise-like quality. The Spectral
Flatness Measure (SFM) is used to compute the tonality of each frame. It is
defined as the ratio of the geometric mean to the arithmetic mean of the power
spectrum P, as per equation (3) to (5).

Figure 4: Tonality measurement and it’s spectrogram for genre disco Figure 5:
Tonality measurement and it’s spectrogram for genre reggae

$\displaystyle P(k)={\rm Re}^{2}[X(k)]+{\rm Im}^{2}[X(k)]$ (3) $\displaystyle
SFMdB=10\log_{10}{\frac{GM\\{P(k)\\}}{AM\\{P(k)\\}}}$ (4) $\displaystyle
Tonality=\min({\frac{SFMdB}{SFMdB_{\max}}},1)$ (5)

####  III-B4 Spectral Centroid

The spectral centroid can be defined as the mean of the distribution of
frequency components for a given frame of the signal. This mean can be
calculated using either the linear frequency or the Bark-scale as parameters.
The weights for each parameter (magnitude of FFT components) are applied
according to Eq. (6)

$\displaystyle
SC=\frac{\sum\limits_{k=0}^{N-1}kX^{2}(k)}{\sum\limits_{k=0}^{N-1}X^{2}(k)}$
(6) $\displaystyle
SC_{b}=\frac{\sum\limits_{j=0}^{N-1}b_{j}(b_{j}-b_{j-1})X^{2}(j)}{\sum\limits_{j=0}^{N-1}(b_{j}-b_{j-1})X^{2}(j)}$
(7)

Figure 6: Spectograms of spectral centroids of two different genres of music
samples

The signal’s sound is affected by its spectral centroid. A higher spectral
centroid indicates a brighter, happier sound, while a lower spectral centroid
indicates a duller, gloomier sound. This is evident in Figure 5. The spectral
centroid, computed over the Bark scale, is a psycho-acoustically adopted a
measure that indicates this Eq. (7).

####  III-B5 Harmonicity

Harmonicity features are a set of characteristics used to analyze the periodic
properties of a signal. These features are based on two primary measures: the
harmonicity ratio and the fundamental frequency. The harmonicity ratio is a
metric that reflects how regularly the signal oscillates, while the
fundamental frequency is the frequency that gives the most coherent
explanation of the signal’s spectrum. The fundamental frequency is computed
using Goldstein’s algorithm [12], which utilizes a likelihood approximation
method to obtain the fundamental frequency.

####  III-B6 Mel-Frequency Cepstral Coefficients

The MFCC, or Mel-Frequency Cepstral Coefficients [13], is a method for
representing the shape of a spectrum using a limited number of coefficients.
This method is based on the cepstrum, which is the Fourier transform of the
logarithm of the spectrum. However, the MFCC uses a variation of the cepstrum
that is calculated on Mel-frequency bands rather than the traditional Fourier
spectrum. This variation, known as the Mel-cepstrum, is particularly effective
at capturing the characteristics of the mid-frequency range of a signal. The
calculation of the Mel-cepstrum is described by equation (8).

$f_{mel}=2595\log_{10}\left(1+\frac{f}{700}\right)$ (8)

####  III-B7 Time Domain Zero Crossings

Zero-crossings in the time domain represents the noisiness of the signal. It
is calculated by using the sign function: 0 for negative arguments while a
positive argument is given for 1 in the signal. Let’s take a signal x[n] in
the time domain. The time domain zero crossings are calculated for the frame t
as per Eq. (9).

$TDZC_{t}=\frac{1}{2}\sum_{n=1}^{M}|sign[x[n]]-sign[x[n-1]]|$ (9)

###  III-C Preprocessing of Dataset

The GTZAN dataset [4] has been used for the work. The GTZAN dataset is a
public domain dataset that consists of 1000 music signals. The music signals
are of 30 seconds in duration. The music signals are divided into 10 genres.
The genres are blues, classical, country, disco, hip-hop, jazz, metal, pop,
reggae, and rock. The GTZAN dataset is divided into a training set and a test
set. The training set consists of 800 music signals and the test set consists
of 200 music signals.

###  III-D Classification Algorithm

The classification of music genres is a challenging task due to the inherent
variability and subjectivity of music. In this study, we proposed a machine-
learning algorithm for music genre classification using the GTZAN dataset. The
algorithm is implemented using the Python programming language and several
libraries such as numpy, pandas, os, librosa, sklearn, and keras.

Algorithm 1 Music Signal Feature Classification

1:procedure Classify($x(t)$)

2: $y(t)\leftarrow$ Preprocess($x(t)$)

3: $f_{1},f_{2},\dots,f_{n}\leftarrow$ ExtractFeatures($y(t)$)

4: $c\leftarrow$ TrainClassifier($f_{1},f_{2},\dots,f_{n}$)

5: return $c(f_{1},f_{2},\dots,f_{n})$

6:end procedure

7:procedure Preprocess($x(t)$)

8: $y(t)\leftarrow$ Downsample($x(t)$)

9: $y(t)\leftarrow$ RemoveNoise($y(t)$)

10: return $y(t)$

11:end procedure

12:procedure ExtractFeatures($y(t)$)

13: $f_{1}\leftarrow$ ComputeMFCC($y(t)$)

14: $f_{2}\leftarrow$ ComputeChroma($y(t)$)

15: $f_{3}\leftarrow$ ComputeSpectralContrast($y(t)$)

16: $\dots$

17: return $f_{1},f_{2},\dots,f_{n}$

18:end procedure

19:procedure TrainClassifier($f_{1},f_{2},\dots,f_{n}$)

20: $c\leftarrow$ SVM($f_{1},f_{2},\dots,f_{n}$)

21: return $c$

22:end procedure

The dataset consists of 1000 audio files of 10 different music genres, with
100 samples per genre. The dataset was preprocessed by extracting the
filenames of the audio files, extracting the labels of the audio files, and
encoding the labels using the LabelEncoder. The labels were then converted to
a categorical format. The features of the audio files were extracted using the
librosa library, which is a library for music and audio processing in Python.
The Mel-Frequency Cepstral Coefficients (MFCCs) were used as the feature
representation of the audio files. The MFCCs were extracted from the audio
data and the mean of the MFCCs was taken across the time axis. The feature
data was then converted to a numpy array. The classification model was built
using the Keras library, which is a high-level neural networks API, written in
Python and capable of running on top of TensorFlow. The model was implemented
as a sequential model with two dense layers. The first dense layer had 256
neurons and the activation function used was ReLU. A 0.5 dropout rate was used
for reducing overfitting. The second dense layer had 9 neurons, corresponding
to the number of genres in the dataset, and the activation function used was
softmax.

Figure 7: Block diagram of the classification system  Figure 8: Deep-Learning
Model of the Classification System

The model was compiled using the categorical cross-entropy loss function and
the Adam optimizer. The model was trained using the feature data and labels,
with a batch size of 40 and 20 epochs. The validation split was set to 10% to
evaluate the performance of the model on unseen data during training. The
proposed algorithm achieved an overall accuracy of 81% on the test data,
demonstrating its effectiveness in classifying music genres. The algorithm can
be further improved by using different feature representations, and by using
more advanced neural network architectures such as convolutional neural
networks.

###  III-E Results and Discussion

The trained model was used to predict the test data and the classification
report was generated using the metrics library. The classification report
provides the precision, recall, f1-score, and support for each genre, which
are useful in evaluating the performance of the model. The diagonal elements
of the confusion matrix are highlighted in the heatmap, where the darker color
represents the higher count of correctly classified observations. The off-
diagonal elements represent the misclassification, whereas the lighter color
represents the lower count of misclassification.

TABLE I: Classification report Music Genre | precision | recall | f1-score | support
---|---|---|---|---
blues | 0.73 | 0.90 | 0.80 | 100
classical | 0.96 | 0.99 | 0.98 | 100
country | 0.75 | 0.94 | 0.83 | 100
disco | 0.63 | 0.85 | 0.73 | 100
hiphop | 0.83 | 0.79 | 0.81 | 100
metal | 0.94 | 0.96 | 0.95 | 100
pop | 0.82 | 0.81 | 0.81 | 100
reggae | 1.00 | 0.07 | 0.13 | 100
rock | 0.98 | 0.06 | 0.87 | 100
Accuracy | 0.80
macro avg | 0.83 | 0.80 | 0.77 | 900
weighted avg | 0.83 | 0.80 | 0.77 | 900
TABLE II: Comparison of the proposed model with available models Sl no. | Model Name | Accuracy |  |
---|---|---|---|---
1 | Gaussian Model | 61% |  |
2 | Logistic Regression Model | 75% |  |
3 | CRNN Model | 53.5% |  |
4 | CNN-RNN Model | 56.4% |  |
5 | Simple Artificial Neural Network | 64.06% |  |
6 | Proposed Model | 80% |  |
Figure 9: Genre classification confusion matrix

It is also important to note that the accuracy of the model can be calculated
using the formula (correctly classified observations / total observations) and
it can be computed using the diagonal elements of the matrix.

##  IV Deployment Architecture

Deployment is an essential step in the development of any machine learning
system. In this section, we propose a deployment architecture for a music
genre classification system that utilizes the cloud services provided by
Amazon Web Services (AWS). The proposed architecture is designed to be
scalable, durable, and easy to access for users. The first component of the
proposed architecture is Amazon S3 [14]. S3 is a fully managed object storage
service that provides scalable and durable storage for audio files and
metadata. This allows for easy management and retrieval of the data needed for
training and inference. Amazon SageMaker [15] is the second component of the
proposed architecture. SageMaker provides a fully managed platform for
building, training, and deploying machine learning models. With SageMaker, we
can train a model for music genre classification using the audio files and
metadata stored in S3. Once the model is trained, it will be deployed to a
SageMaker endpoint for inference. The endpoint can be accessed via an API,
allowing users to submit audio files for classification. Amazon API Gateway
[16] is used to create a RESTful API for the SageMaker endpoint, providing a
convenient way for users to access the classification service. The
classification results and metadata will be stored in Amazon DynamoDB [17], a
fully managed NoSQL database service. DynamoDB provides high performance and
scalability, making it well-suited for storing large amounts of data generated
by a music genre classification system.

Figure 10: Block diagram of the proposed Deployment Architecture for
Mobile/Web apps

To enable fast and powerful search capabilities for the classification results
and metadata, an Elasticsearch index will be created using Amazon
Elasticsearch (OpenSearch [18]) Service. Elasticsearch is a popular search
engine that is well-suited for handling large amounts of data. Finally, Amazon
CloudFront [19] will be used to distribute the classification results and
metadata to users. CloudFront is a content delivery network (CDN) that ensures
low latency and high availability of the results, making it easy for users to
access the classification results from anywhere in the world. In conclusion,
the proposed architecture is designed to provide a scalable, durable, and
easy-to-access music genre classification system using the cloud services
provided by AWS. The architecture includes various services like S3,
SageMaker, API Gateway, DynamoDB, Elasticsearch, and CloudFront. These
services together enable the system to handle a large amount of data, train
and deploy models effectively and provide fast and accurate

##  V Conclusion

The proposed DSP-based music genre classification system was found to be
effective in classifying various types of music. The system was able to
correctly classify different types of music with an accuracy of 80%. The
proposed system can be used to classify different types of music in a real-
time scenario and also when the music was played at different speeds. The
proposed system can also be used to automatically generate playlists for users
based on their music preferences and it will help researchers to better
understand the relationship between music and human emotions.

## References

  * [1] A. S. Bhat, V. Amith, N. S. Prasad, and D. M. Mohan, “An efficient classification algorithm for music mood detection in western and hindi music using audio feature extraction,” in 2014 Fifth International Conference on Signal and Image Processing, IEEE, Jan. 2014.
  * [2] J. Kim, S. Lee, and W. Yoo, “Implementation and analysis of mood-based music recommendation system,” in 2013 15th International Conference on Advanced Communications Technology (ICACT), pp. 740–743, 2013.
  * [3] J. K. Patel and E. Gopi, “Musical notes identification using digital signal processing,” Procedia Computer Science, vol. 57, pp. 876–884, 2015.
  * [4] G. Tzanetakis and P. Cook, “Musical genre classification of audio signals,” in Proceedings of the International Conference on Acoustics, Speech, and Signal Processing (ICASSP), pp. 41–44, 2002.
  * [5] E. E. P. Myint and M. Pwint, “An approach for mulit-label music mood classification,” in 2010 2nd International Conference on Signal Processing Systems, IEEE, July 2010.
  * [6] M. Shah, G. Wichern, A. Spanias, and H. Thornburg, “Audio content-based feature extraction algorithms using j-DSP for arts, media and engineering courses,” in 2010 IEEE Frontiers in Education Conference (FIE), IEEE, Oct. 2010.
  * [7] S. Dara and P. Tumma, “Feature extraction by using deep learning: A survey,” in 2018 Second International Conference on Electronics, Communication and Aerospace Technology (ICECA), IEEE, Mar. 2018.
  * [8] J. Zhang, “Music feature extraction and classification algorithm based on deep learning,” Scientific Programming, vol. 2021, pp. 1–9, May 2021.
  * [9] K. Zhang, “Music style classification algorithm based on music feature extraction and deep neural network,” Wireless Communications and Mobile Computing, vol. 2021, pp. 1–7, Sept. 2021.
  * [10] S. B. Alex, L. Mary, and B. P. Babu, “Attention and feature selection for automatic speech emotion recognition using utterance and syllable-level prosodic features,” Circuits, Systems, and Signal Processing, vol. 39, pp. 5681–5709, May 2020.
  * [11] M. D. Pawar and R. D. Kokate, “Convolution neural network based automatic speech emotion recognition using mel-frequency cepstrum coefficients,” Multimedia Tools and Applications, vol. 80, pp. 15563–15587, Feb. 2021.
  * [12] S. K. Yoon and S. kyu Kim, “Modified goldstein algorithm using boundary information in phase unwrapping,” in Advances in Imaging, OSA, 2009.
  * [13] M. G. de Pinto, M. Polignano, P. Lops, and G. Semeraro, “Emotions understanding model from spoken language using deep neural networks and mel-frequency cepstral coefficients,” in 2020 IEEE Conference on Evolving and Adaptive Intelligent Systems (EAIS), IEEE, May 2020.
  * [14] M. R. Palankar, A. Iamnitchi, M. Ripeanu, and S. Garfinkel, “Amazon s3 for science grids: A viable solution?,” in Proceedings of the 2008 International Workshop on Data-Aware Distributed Computing, DADC ’08, (New York, NY, USA), p. 55–64, Association for Computing Machinery, 2008.
  * [15] E. Liberty, Z. Karnin, B. Xiang, L. Rouesnel, B. Coskun, R. Nallapati, J. Delgado, A. Sadoughi, Y. Astashonok, P. Das, C. Balioglu, S. Chakravarty, M. Jha, P. Gautier, D. Arpin, T. Januschowski, V. Flunkert, Y. Wang, J. Gasthaus, L. Stella, S. Rangapuram, D. Salinas, S. Schelter, and A. Smola, “Elastic machine learning algorithms in amazon sagemaker,” in Proceedings of the 2020 ACM SIGMOD International Conference on Management of Data, SIGMOD ’20, (New York, NY, USA), p. 731–737, Association for Computing Machinery, 2020.
  * [16] A. Baird, S. Buliani, V. Nagrani, and A. Nair, “Aws serverless multi-tier architectures; using amazon api gateway and aws lambda,” Amazon Web Services Inc, 2015.
  * [17] S. Sivasubramanian, “Amazon dynamodb: A seamlessly scalable non-relational database service,” in Proceedings of the 2012 ACM SIGMOD International Conference on Management of Data, SIGMOD ’12, (New York, NY, USA), p. 729–730, Association for Computing Machinery, 2012.
  * [18] V.-A. Zamfir, M. Carabas, C. Carabas, and N. Tapus, “Systems monitoring and big data analysis using the elasticsearch system,” in 2019 22nd International Conference on Control Systems and Computer Science (CSCS), pp. 188–193, 2019.
  * [19] A. Toshniwal, K. S. Rathore, A. Dubey, P. Dhasal, and R. Maheshwari, “Media streaming in cloud with special reference to amazon web services: A comprehensive review,” in 2020 4th International Conference on Intelligent Computing and Control Systems (ICICCS), pp. 368–372, 2020.


# Impact of periodic vaccination in SEIRS seasonal model

Enrique C. Gabrick<EMAIL_ADDRESS>
<EMAIL_ADDRESS>Potsdam Institute for Climate Impact Research,
Telegrafenberg A31, 14473 Potsdam, Germany.  Department of Physics, Humboldt
University Berlin, Newtonstraße 15, 12489 Berlin, Germany.  Graduate Program
in Science, State University of Ponta Grossa, 84030-900, Ponta Grossa, PR,
Brazil.  Eduardo L. Brugnago  Institute of Physics, University of São Paulo,
05508-090, São Paulo, SP, Brazil.  Silvio L.T. de Souza  Federal University of
São João del-Rei, Campus Centro-Oeste, 35501-296, Divinópolis, MG, Brazil.
Kelly C. Iarosz  Graduate Program in Science, State University of Ponta
Grossa, 84030-900, Ponta Grossa, PR, Brazil.  University Center UNIFATEB,
84266-010, Telêmaco Borba, PR, Brazil.  José D. Szezech Jr  Graduate Program
in Science, State University of Ponta Grossa, 84030-900, Ponta Grossa, PR,
Brazil.  Department of Mathematics and Statistics, State University of Ponta
Grossa, 84030-900, Ponta Grossa, PR, Brazil.  Ricardo L. Viana  Institute of
Physics, University of São Paulo, 05508-090, São Paulo, SP, Brazil.
Department of Physics, Federal University of Paraná, 81531-980, Curitiba, PR,
Brazil.  Iberê L. Caldas  Institute of Physics, University of São Paulo,
05508-090, São Paulo, SP, Brazil.  Antonio M. Batista  Graduate Program in
Science, State University of Ponta Grossa, 84030-900, Ponta Grossa, PR,
Brazil.  Institute of Physics, University of São Paulo, 05508-090, São Paulo,
SP, Brazil.  Department of Mathematics and Statistics, State University of
Ponta Grossa, 84030-900, Ponta Grossa, PR, Brazil.  Jürgen Kurths  Potsdam
Institute for Climate Impact Research, Telegrafenberg A31, 14473 Potsdam,
Germany.  Department of Physics, Humboldt University Berlin, Newtonstraße 15,
12489 Berlin, Germany.

###### Abstract

We study three different strategies of vaccination in a SEIRS
(Susceptible–Exposed–Infected–Recovered–Susceptible) seasonal forced model,
which are: ($i$) continuous vaccination; ($ii$) periodic short time localized
vaccination and ($iii$) periodic pulsed width campaign. Considering the first
strategy, we obtain an expression for the basic reproduction number and infer
a minimum vaccination rate necessary to ensure the stability of the disease-
free equilibrium (DFE) solution. In the second strategy, the short duration
pulses are added to a constant baseline vaccination rate. The pulse is applied
according to the seasonal forcing phases. The best outcome is obtained by
locating the intensive immunization at inflection of the transmissivity curve.
There, a vaccination rate of $44.4\%$ of susceptible individuals is enough to
ensure DFE. For the third vaccination proposal, additionally to the amplitude,
the pulses have a prolonged time width. We obtain a non-linear relationship
between vaccination rates and the duration of the campaign. Our simulations
show that the baseline rates, as well as the pulse duration, can substantially
improve the vaccination campaign effectiveness. These findings are in
agreement with our analytical expression. We show a relationship between the
vaccination parameters and the accumulated number of infected individuals,
over the years and show the relevance of the immunisation campaign annual
reaching for controlling the infection spreading. Regarding the dynamical
behaviour of the model, our simulations shows that chaotic and periodic
solutions, as well as bi-stable regions, depend on the vaccination parameters
range.

> The spread of infectious diseases is a challenge to world public health.
> Many efforts are dedicated to mitigating the impacts of the spreading of
> diseases. In this context, mathematical model is a powerful tool to
> simulate, forecast and study efficient control measures for human and
> wildlife diseases. Although there are many types of disease, some of them
> have common characteristics, for instance to repeat the peak of infectious
> in certain times. These diseases are called seasonal, e.g., measles, mumps
> and smallpox. The reasons for the seasonality can be varied, such as climate
> and social behaviours. Due to the seasonal characteristic of these diseases,
> they recur in the populations and control measures are needed to be
> implemented in order to eradicated them. One the most successful control
> measure is a vaccination campaign, that can be continuous or periodic.
> Continuous vaccination is a campaign in which a constant quantity of
> vaccines is daily available to the population throughout the year, while the
> periodic strategy considers an intense immunization campaign, during a
> fraction of the year. In this work, we study the impacts of vaccination in a
> SEIRS model with seasonal forcing. We consider newborns vaccine and
> susceptible vaccination, administrated in three different ways: ($i$)
> continuous vaccination; ($ii$) periodic short time localized vaccination and
> ($iii$) periodic pulsed width campaign. For the periodic strategies, we
> consider a baseline rate. In the constant vaccine context, we obtain an
> analytical expression for the vaccination rate in which the (DFE) is
> guaranteed. For the parameters considered, the annually vaccinated value is
> in agreement with numerical simulations for all the vaccine strategies.
> Furthermore, considering parameters for which bi-stability exist in the case
> without vaccine, our results show that it persists depending on the
> vaccination rate.

##  I Introduction

Infectious diseases spreading is a highly important problem in public health
Xia2023 . The spread of infectious diseases affects the humanity throughout
whole history Glatter2021 , e.g., Black Death during the fourteenth century
Glatter2021 , Spanish flu in 1918 Tumpey2005 , COVID-19 Machein2020 , Dengue
Fever Aguiar2011 , HIV Dalal2008 , etc. Meyers2005 ; Mao2010 ; Scarpino2019 .
Some of these illness present seasonal patterns of incidence Altizer2006 ,
namely seasonal infections diseases Buonomo2018 . The mechanisms of
seasonality are varied, such as weather Hoshen2004 , school holidays
Finkenstadt2000  and others Grassly2006 . Examples of seasonal infectious
diseases are mumps Azimaqin2022 , measles Xia2023  and dengue fever Aguiar2009
. An introduction to corresponding models is found in Refs. Altizer2006 ;
Buonomo2018 .

Mathematical models are a powerful tool to understand, forecast and study
diseases spread control measures Mugnaine2022 . In epidemiology, the standard
models are compartmental Batista2021 , where the host population is
compartmentalized according to the stages of the infection evolution and the
possible states considered. The individuals in the host population are taken
from susceptible to possible intermediate stages, as infectious and recovered.
Depending on the disease to be modelled, when an infection cycle is completed,
the individuals can acquire permanent immunity Grenfell1995 , or become
susceptible again. These models are easily adapted to study different
diseases, as seasonal infectious ones Keeling2008 , which can be modelled by
the inclusion of a nonlinear time-dependent term in the transmissivity, e.g.,
a sinusoidal forcing Bilal2016  or a square wave function Tanaka2013 . Due to
this non-linearity, the resulting dynamics can become very intricate, even
exhibiting chaos Olsen1990  or bi-stability Gabrick2023 .

From a modelling perspective, seasonal models have been used since 1928
Buonomo2018 . These models reproduce with great accuracy the dynamics observed
in diseases, such as measles Olsen1990 , chickenpox, mumps London1973  and
others Altizer2006 . In a formulation of the SIRS
(Susceptible–Infected–Recovered–Susceptible) model with seasonal contact rate,
Greenhalgh and Moneim Greenhalgh2003  showed the existence of an unique DFE
solution, which is globally asymptotically stable when the basic reproduction
number is less than one ($\mathcal{R}_{0}<1$), where $\mathcal{R}_{0}$ is a
measure of the reproductive potential for a given disease. In a population,
where everyone is initially susceptible, the infection can remain only if
$\mathcal{R}_{0}>1$ Keeling2008 . Furthermore, Greenhalgh and Moneim
considered four childhood infectious diseases (measles, chickenpox, mumps, and
rubella) and showing non-trivial solutions (chaotic dynamics). These results
are found by analysing the bifurcation diagram, where some ranges exhibit
chaotic attractors. Also in a SIR (Susceptible–Infected–Recovered) framework,
de Carvalho and Rodrigues Carvalho2022  implemented a multi-parameter
periodically forced term leading to strange attractor solutions. In their
formulation, the DFE is not preserved when $\mathcal{R}_{0}<1$. Metcalf et al.
Metcalf2009  studied the seasonal variation of six childhood infections
(measles, pertussis, mumps, diphtheria, varicella and scarlet fever) from data
of Copenhagen in the pre-vaccination era. Their results showed that the
transmission disease decreases for some infections at school holidays.

Previous works reported the capacity of compartmental models to reproduce
chaotic dynamics Keeling2001  and the accuracy in simulating real data
Machein2020 ; Cooper2020 . However, one of the most important advantages in
using compartmental models is the facility to implement important
characteristics to simulate control measures, such as social distancing
Ansari2021 ; Mello2020 , restrictive measures Mugnaine2022 ; Brugnago2020 ;
Souza2021 , quarantine Balsa2021 ; Zou2022 , vaccination Etxeberria-
Etxaniz2020 ; DelaSen2021 , etc. Despite there are some forms of control, one
of the most effective is vaccination Voysey2021 ; Thompson2022 ; Zou2022 .

Vaccination campaign in infectious seasonal diseases shows a significant
infected number decrease Grassly2006 ; Molinari2007 ; Keeling2002 . Gao et al.
Gao2011  studied the vaccination of newborns combined with pulsed vaccine in a
SIRS seasonal forced model, in a modelling approach. One of the results
obtained in this research is the DFE when $\mathcal{R}_{0}<1$. The rotavirus
vaccination was studied by Atchison et al. Atchison2010  using a modified SIR
model to fit the data from England and Wales. From these simulations, they
reported that vaccination reduces rotavirus diseases transmission by $61\%$.
Considering a pulsed vaccination strategy in a SIR model, Shulgin et al.
Shulgin1998  showed the possibility of disease eradication with relatively low
vaccination rates. To get the disease eradication, they explored some
conditions, as vaccine proportion and periodicity. The effects of two
vaccination doses in a SEIR (Susceptible-Exposed-Infected-Recovered) epidemic
model was studied by Gabrick et al. Gabrick2022 . In this work, they
considered three vaccination strategies: unlimited doses applied continuously
in the susceptible population; limited doses supply applied periodically in
the susceptible population; and limited doses applied in a periodic strategy
in all host population. Their results showed that the vaccine campaign is more
efficient when applied only in susceptible individuals, i.e., the population
is previously tested. Considering a Kot-type function as seasonality, Duarte
Duarte2021  analysed the control of infectious disease with vaccination
strategies including a perturbation term. Seasonal contact and optimal
vaccination strategy were studied by Wang Wang2019  in a SEIR model. The
persistence or extinction of a seasonal disease with reinfection possibility
was discussed by Bai and Zhou Bai2012 . They explored the conditions for the
extinction of the diseases in the situations where $\mathcal{R}_{0}<1$ and,
also $\mathcal{R}_{0}>1$. Their results showed that $\mathcal{R}_{0}=1$ is the
threshold for the disease extinction. However, from simulations, a policy only
based on $\mathcal{R}_{0}$ can overestimate the infections risks and the
infected number due the presence of seasonality. Periodic strategy vaccination
in a SEIRS model was discussed by Moneim and Greenhalgh Moneim2005 . The
vaccine periodicity was described as integer multiples of the contact rate
period. They reported that a key parameter to understand the vaccine influence
is the $\mathcal{R}_{0}$. Also, they proved that $\mathcal{R}_{0}<1$ is
associated with the DFE and is globally asymptotically stable. Other works
that address the vaccination of seasonal diseases can be found in Refs. Ho2019
; Metcalf2012 ; Nokes1995 ; Agur1993  and in the references therein.

In this work, we study the effects of vaccination in a SEIRS seasonal forced
model Bai2012 ; Gabrick2023 . We substantially extended the works of Moneim
and Greenhalgh Moneim2005  and Gabrick et al. Gabrick2023 . The first one
explored the effects of periodic vaccination and conjectured one
$\mathcal{R}_{0}$ as a function of vaccine parameters. However, the authors
did not explore the effects of constant vaccine and periodic pulses with
different width. In the second one, the authors obtained a parameter range in
which the numerical solutions exhibit bi-stability and tipping point phenomena
can be explored. However, they did not taken into account the vaccine
influences in the dynamical behaviour. In this way, our work extend the Ref.
Moneim2005  by implementing three different strategies of vaccination in
susceptible and newborns individuals, namely ($i$) constant, ($ii$) pulsed and
($iii$) pulsed width vaccination strategies. We choose these three strategies
because: ($i$) it can be used to model mass vaccination programs
Grabenstein2006 ; ($ii$) it can be considered to model a scenario in which a
certain amount of vaccination (baseline) is available during all the time, but
periodically an intense campaign are employed Alkhamis2022 ; and strategy
($iii$) modified ($ii$) by extending the duration of the intensive
immunization campaign to days, weeks or months, i.e., it saves vaccination
efforts. Then, considering the parameters found in Ref. Gabrick2023  for bi-
stability, we explore, as novelty, the dynamical aspects of the SEIRS seasonal
in the presence of vaccine, according to the three strategies mentioned. Our
simulations show that for the ($i$) case the vaccine in susceptible is more
significant and rates $\geq 44\%$ are able to extinct the infection in the
host population. In the ($ii$) scenario we found a linear relationship between
the baseline and pulsed immunization rate to result in the infection
eradication. Also, acting the pulses at the inflection point of the
seasonality function, the campaign is slightly more effective. To the strategy
($iii$), we discover a non-linear relationship between the immunization
campaign parameters that leads to the illness extinction.

The current work is organised as follows: in Sec. II we present the model and
its equilibrium solutions, from which we obtain the $\mathcal{R}_{0}$ as well
as a minimum rate to eradicate the disease. In Sec. III we discuss the effects
of constant vaccination in the model. The pulsed vaccination strategy is
discussed in Sec. IV and the implementation of pulsed width in the vaccine
campaign is present in Sec. V. Finally, our conclusions are drawn in Sec. VI.

##  II Model

SEIRS is a compartmental epidemiological model that describes the spread of a
given infectious disease in a homogeneously mixed host population, in which
individuals are computed into one of four compartments, namely susceptible in
$S$, exposed and not yet contagious in $E$, infectious in $I$ and recovered in
$R$ Keeling2008 . The population size ($N=S+E+I+R$) is time-dependent when the
natural death rate ($\mu$) is not equal to the birth rate ($b$). The
transition flow between the compartments is schematically represented in Fig.
1. Susceptible individuals become infected through interaction with contagious
agents at a rate $\beta{I/N}$, where $\beta$ is the transmitting infection
rate per interaction. Once a portion of $S$ is infected, it evolves to $E$,
remaining in this compartment by an average time given by $1/\alpha$ (latent
period). The parameter $\alpha$ is the rate at which exposed individuals
become infectious, migrating to compartment $I$. Individuals remain in $I$ by
an average time $1/\gamma$ (infectious interval), after that they occupy the
$R$ compartment, the parameter $\gamma$ is called recovery rate. In this model
there is no permanent immunity, then, after an average interval $1/\delta$ the
individuals in $R$ return to $S$ Batista2021 . In which $\delta$ is the
immunity loss rate. Furthermore, we consider vaccination of newborns and
susceptible Moneim2005 . A fraction $p\in[0,1]$ of newborns are directly
immunized and a portion of $S$ is vaccinated at a rate $v$ with effectiveness
$\lambda\in[0,1]$. The vaccine effect is to give immunity to individuals,
transferring, in this model, newborn and susceptible ones to the recovered
class. Given the nature of the model, all parameters are non-negative real
numbers.

Figure 1: Schematic representation of SEIRS model with vaccine inclusion. The
compartments (coloured boxes), indexed by variables of the system (1),
contains fractions of the population at each stage of infection spread. Arrows
indicate the flow between them, increase in newborns and loss of population
due to death; accompany the respective transition rates.

The SEIRS model with vaccination, as described above, is given by the
following system of four coupled ordinary differential equations Moneim2005 :

$\begin{split}\frac{dS}{dt}&=b(1-p)N+\delta R-\left(\beta\frac{I}{N}+\lambda
v+\mu\right)S,\\\ \frac{dE}{dt}&=\beta\frac{SI}{N}-(\alpha+\mu)E,\\\
\frac{dI}{dt}&=\alpha E-(\gamma+\mu)I,\\\ \frac{dR}{dt}&=bpN+\lambda vS+\gamma
I-(\delta+\mu)R.\end{split}$ (1)

Other effects also can be considered in this model and incorporated into Eqs.
(1), such as reaction-diffusion Song2019  and delay Cooke1996 , but this is
beyond the present work and will be studied next. Seasonality is included in
the model by replacing $\beta$ by a periodic function:

$\beta(t)=\beta_{0}\left[1+\beta_{1}{\rm cos}(\omega t)\right],$ (2)

where $\beta_{0}$ is the average contagion rate, the seasonality degree is
$\beta_{1}\in[0,1]$, and $\omega$ is its frequency Olsen1990 .

From Eqs. (1), the sum of variables followed by a simple manipulation,
provides us the growth rate of the host population $dN/dt=(b-\mu)N$.
Considering this result, the Eqs. (1) can be rewritten in a normalised form,
without loss of generality, by taking Bjornstad2018

$S=Ns,~{}E=Ne,~{}I=Ni,~{}\text{and}~{}R=Nr,$

which gives the respective transformations for the time derivatives of the
variables:

$\begin{split}\frac{ds}{dt}&=\frac{1}{N}\left(\frac{dS}{dt}-(b-\mu)S\right),\\\\[1.5pt]
\frac{de}{dt}&=\frac{1}{N}\left(\frac{dE}{dt}-(b-\mu)E\right),\\\\[1.5pt]
\frac{di}{dt}&=\frac{1}{N}\left(\frac{dI}{dt}-(b-\mu)I\right),\\\\[1.5pt]
\frac{dr}{dt}&=\frac{1}{N}\left(\frac{dR}{dt}-(b-\mu)R\right).\end{split}$ (3)

This change of variables leads to equations of the same form that would be
obtained by assuming $b=\mu$. We get

$\begin{split}\frac{ds}{dt}&=b(1-p)+\delta r-\left[\beta(t)i+\lambda
v+b\right]s,\\\ \frac{de}{dt}&=\beta(t)is-(\alpha+b)e,\\\
\frac{di}{dt}&=\alpha e-(\gamma+b)i,\\\ \frac{dr}{dt}&=bp+\lambda vs+\gamma
i-(\delta+b)r,\end{split}$ (4)

with the constrain $s+e+i+r=1$. Then, $r$ can be determined in terms of the
other three variables, such that we replace it in the first of Eqs. (4),
reducing the model to the following three-equation system:

$\displaystyle\frac{ds}{dt}$
$\displaystyle=b(1-p)+\delta-\left[\beta(t)i+\lambda
v+\delta+b\right]s-\delta(e+i),$ $\displaystyle\frac{de}{dt}$
$\displaystyle=\beta(t)is-(\alpha+b)e,$ (5) $\displaystyle\frac{di}{dt}$
$\displaystyle=\alpha e-(\gamma+b)i.$

###  II.1 Disease-free equilibrium

In general, Eqs. (5) are solved numerically. However, for the particular case
when $\beta=\beta_{0}$ is constant, there are two fixed point solutions, which
are: DFE and endemic solution Guillen2017 . The first one is characterised by
the disappearance of the disease in the host population ($i^{*}=0$); in the
other, the infection remains ($i^{*}>0$). Firstly, we investigate the DFE
solution given by

$(s^{*},e^{*},i^{*})=\left(\frac{b(1-p)+\delta}{\lambda
v+\delta+b},0,0\right).$ (6)

This fixed point is stable if all eigenvalues of the Jacobian matrix
($\mathbb{J}$) of the system, computed in DFE, have negative real part
TamasTel . Given the matrix

$\mathbb{J}_{\rm DFE}=-\left[\begin{array}[]{ccc}(\lambda
v+\delta+b)&\delta&(\beta_{0}s^{*}+\delta)\\\\[4.0pt]
0&(\alpha+b)&-\beta_{0}s^{*}\\\\[4.0pt]
0&-\alpha&(\gamma+b)\end{array}\right],$ (7)

its eigenvalues are

$\displaystyle\xi_{1}$ $\displaystyle=-(\lambda v+\delta+b),$ (8)
$\displaystyle\xi_{2}$
$\displaystyle=\frac{-(\alpha+\gamma+2b)-\sqrt{(\alpha-\gamma)^{2}+4\alpha\beta_{0}s^{*}}}{2},$
(9) $\displaystyle\xi_{3}$
$\displaystyle=\frac{-(\alpha+\gamma+2b)+\sqrt{(\alpha-\gamma)^{2}+4\alpha\beta_{0}s^{*}}}{2}.$
(10)

All of them are real numbers, once the constants are positive and $\xi_{1,2}$
are always negative, then the DFE is stable when $\xi_{3}<0$. From this
inequality, we obtain the relation for the stability of this fixed point:

$\mathcal{R}_{0}=\frac{\alpha\beta_{0}s^{*}}{(\alpha+b)(\gamma+b)}<1.$ (11)

Therefore, if $\mathcal{R}_{0}<1$ the DFE point is stable, i.e. the disease
will die out. On the other hand, if $\mathcal{R}_{0}>1$ the disease will be
succeed in infect the population AndersonMay . Thus, we establish a stability
condition for the DFE relative to the basic reproduction number. The quantity
$\mathcal{R}_{0}$ informs about the evolution of the spread. Note that, if
$p=\lambda v=0$, then $s^{*}=1$ and we recover the $\mathcal{R}_{0}$ for the
model without vaccine Gabrick2023 .

Note that Eq. (11) is for the autonomous case, with $\beta$ and $v$ is
constant. However, for periodic diseases without latent period Allen2010  the
$\mathcal{R}_{0}$ expression can be extend replacing $\beta_{0}$ by
$\langle{\beta(t)}\rangle$:

$\langle{\beta(t)}\rangle=\frac{1}{T_{\rm f}-T_{0}}\int_{T_{0}}^{T_{\rm
f}}\beta(t)dt,$ (12)

where $T_{\rm f}-T_{0}$ corresponds to the seasonality period. In the case
studied in this work, we can analyze both the average
$\langle\mathcal{R}\rangle$ of the time-dependent reproduction number and its
maximum $\mathcal{R}_{+}$ within a seasonal cycle Moneim2005 , being
$\langle\mathcal{R}\rangle=\mathcal{R}_{0}$ and
$\mathcal{R}_{+}=(1+\beta_{1})\mathcal{R}_{0}$. imposing the stable DFE
condition on these quantities, we obtain two quotas:

$\displaystyle\frac{\alpha[b(1-p)+\delta]}{\lambda(\alpha+b)(\gamma+b)}\beta_{0}-\frac{(\delta+b)}{\lambda}$
$\displaystyle<v_{0},$ (13)
$\displaystyle\frac{\alpha[b(1-p)+\delta]}{\lambda(\alpha+b)(\gamma+b)}\beta_{0}(1+\beta_{1})-\frac{(\delta+b)}{\lambda}$
$\displaystyle<v_{+}.$ (14)

The disease extinction is guaranteed when the relationship $v_{+}\leq v$ is
verified. However, this upper limit for the vaccination rate can be greater
than enough, with $v_{0}\leq v$ being sufficient to reach the DFE in our
simulations. Note that this threshold is the same in the autonomous case. As a
numerical example, considering $\lambda=1$, $p=0.25$ $b=0.02$,
$\beta_{0}=270$, $\delta=0.25$, $\alpha=100$ and $\gamma=100$, we obtain a
quota $v_{0}=0.445$. Therefore, for a constant vaccination campaign the
minimum rate is $44.5\%$ of vaccinated susceptible. A value very close to
those numeriacally obtained in Sec. III and annual vaccination rates in Secs.
IV and V.

###  II.2 Endemic equilibrium

Considering the autonomous case, the endemic equilibrium point ${\rm
EE}\left(\widetilde{s},\widetilde{e},\widetilde{i}\right)$ is given in terms
of $\mathcal{R}_{0}$ and $s^{*}$, as

$\displaystyle\widetilde{s}$
$\displaystyle=\frac{(\alpha+b)(\gamma+b)}{\alpha\beta_{0}}=\frac{s^{*}}{\mathcal{R}_{0}},$
(15) $\displaystyle\widetilde{e}$
$\displaystyle=\frac{(\gamma+b)\left[b(1-p)+\delta\right]}{\alpha\beta_{0}s^{*}+\delta(\alpha+\gamma+b)\mathcal{R}_{0}}\left(\mathcal{R}_{0}-1\right),$
(16) $\displaystyle\widetilde{i}$
$\displaystyle=\frac{\alpha\left[b(1-p)+\delta\right]}{\alpha\beta_{0}s^{*}+\delta(\alpha+\gamma+b)\mathcal{R}_{0}}\left(\mathcal{R}_{0}-1\right).$
(17)

In this solution the disease is permanent in the host population. It is worth
to note that for the inverse relation $\widetilde{s}$ proportional to
$\mathcal{R}_{0}^{-1}$, the scenario without vaccination verifies the result
Gabrick2023  $\widetilde{s}=\mathcal{R}_{0}^{-1}$. This fixed point exists
only if $\mathcal{R}_{0}>1$, since $\mathcal{R}_{0}<1$ implies
$\widetilde{e},\widetilde{i}<0$ in Eqs. (16) and (17). It means that, in the
case of stable DFE, there is no endemic fixed point, recovering the
information previously discussed.

As in the DFE case, we analyze the stability of the EE fixed point through
eigenvalues of the Jacobian matrix $\mathbb{J}_{\rm EE}$ computed on it, being

$\mathbb{J}_{\rm
EE}=-\left[\small{\begin{array}[]{ccc}(\beta_{0}\widetilde{i}+\lambda
v+\delta+b)&\delta&(\beta_{0}\widetilde{s}+\delta)\\\\[4.0pt]
-\beta_{0}\widetilde{i}&(\alpha+b)&-\beta_{0}\widetilde{s}\\\\[4.0pt]
0&-\alpha&(\gamma+b)\end{array}}\right].$ (18)

The correspondent characteristic polynomial is given by

$\mathcal{P}(\zeta)=\zeta^{3}+c_{2}\zeta^{2}+c_{1}\zeta+c_{0},$ (19)

where the coefficients are

$\displaystyle c_{0}$ $\displaystyle=(\beta_{0}\widetilde{i}+\lambda
v+\delta+b)\left[(\alpha+b)(\gamma+b)-\alpha\beta_{0}\widetilde{s}\right]+$
$\displaystyle+\beta_{0}\widetilde{i}\left[\alpha\left(\beta_{0}\widetilde{s}+\delta\right)+\delta(\gamma+b)\right],$
(20) $\displaystyle c_{1}$ $\displaystyle=(\beta_{0}\widetilde{i}+\lambda
v+\delta+b)(\alpha+\gamma+2b)+$
$\displaystyle+(\alpha+b)(\gamma+b)+\beta_{0}(\delta\widetilde{i}-\alpha\widetilde{s}),$
(21) $\displaystyle c_{2}$ $\displaystyle=\beta_{0}\widetilde{i}+\lambda
v+\alpha+\gamma+\delta+3b.$ (22)

Having such coefficients and using the Routh–Hurwitz criterion Meinsma1995237
, it is possible to find a parametric relationship for the EE point stability.
In this way, the real part of the all three eigenvalues is less than zero if
and only if

$c_{0},c_{1},c_{2}>0;~{}~{}\text{and}~{}~{}c_{0}<c_{1}c_{2}.$ (23)

Fulfilled the above conditions, once the EE point exists, it is attractive. In
the context of the present work, we verify that the EE point is asymptotically
stable for the parameters used in our simulations, as already described in the
previous subsection: $\lambda=1$, $p=0.25$, $b=0.02$, $\delta=0.25$,
$\alpha=100$, $\gamma=100$, $\beta_{0}=270$ and $v\in[0,v_{0}=0.445)$, where
$\mathcal{R}_{0}>1$.

For the non-autonomous system, where $\beta=\beta(t)$ is the periodic function
given in Eq. (2), the DFE solution is the only one of the fixed-point type.

Figure 2: Trajectories evolving from the initial condition
$P_{0}(s_{0},e_{0},i_{0})=(0.39,0.0039,0.0039)$ (black dot). We consider the
parameters: $v=0.17$, $\lambda=1$, $p=0.25$, $b=0.02$, $\beta_{0}=270$,
$\delta=0.25$, $\alpha=100$ and $\gamma=100$. For $\beta_{1}=0$, the
trajectory (orange line) spirals towards the attractive fixed point ${\rm
EE}(0.371,0.001,0.001)$ (red dot). In the non-autonomous case
($\beta_{1}=0.28$ and $\omega=2\pi$), the system evolves (light-blue line)
around the point ${\rm EE}$ converging to the periodic attractor (blue line).

In other words, the fixed point solution with $i\neq 0$ does not exist when
the transmissivity is an explicit function of time. In this case, endemic
solutions are not fixed points, but consist of orbits that can be periodic or
chaotic. Figure 2 shows the system evolution from the initial condition
$P_{0}(s_{0},e_{0},i_{0})=(0.39,0.0039,0.0039)$ in the two cases: 1)
autonomous (with $\beta_{1}=0$), where the trajectory (orange line) spirals to
the attractive fixed point ${\rm
EE}\left(\widetilde{s},\widetilde{e},\widetilde{i}\right)\approx(0.371,0.001,0.001)$;
2) non-autonomous system (with $\beta_{1}=0.28$), where the trajectory (light-
blue line) evolves during a transient time around EE point and converges to
the periodic attractor (blue line). To obtain these curves, we adopt the
aforementioned parametric configuration and consider a constant vaccination
rate $v=0.17$. In the following sections, we numerically investigate the
impacts of different vaccination protocols on the dynamics of the Eqs. (5)
with seasonality, focusing on the infected population.

##  III Constant vaccination strategy

Figure 3: Numerical solutions of Eqs. (5) for three different values of $p$
and $v$. (a) Limit cycles obtained for $p=v=0$ (red curve) and $p=0.25$, with
$v=0.2$ (blue curve) and $v=0.4$ (green curve). (b) Time series for $i$
corresponding to the limit cycles in the same colours according to the legend.
We consider $b=0.02$, $\beta_{0}=800$, $\beta_{1}=0.20$, $\delta=0.25$,
$\alpha=40$, $\gamma=100$, $\omega=2\pi$, and initial condition
$(s_{0},e_{0},i_{0})=(0.9,0,0.1)$. We discard $10^{5}$ integration steps as
transient.

Continuous vaccination campaign is modelled by a constant vaccine term $v$ in
Eqs. (5), which represents a fraction of the host population being vaccinated
at every time step. This approach can be used for modelled mass vaccination
programs Grabenstein2006 . In this campaign, a certain amount of the
population is vaccinated in a short period of time. Some examples, are the
mass campaign against measles McLean1988 , smallpox Fenner1988 , and COVID-19
Bagcchi2021 . In this case, $p$ is also constant. To numerically integrate the
system of differential equations, we use the fourth-order Runge-Kutta method
Butcher1987  with a fixed step of $10^{-3}$. It is important to mention that,
without spoiling to the analyses, we consider the vaccine effectiveness
$\lambda=1$. Even so, this term can be understood, more generally, as absorbed
by $v$, since these always appear as a product, maintaining the equations
form.

Figure 3 displays three different numerical solutions for Eqs. (5), for
different $p$ and $v$. In these results, we consider $b=0.02$,
$\beta_{0}=800$, $\beta_{1}=0.20$, $\delta=0.25$, $\alpha=40$, $\gamma=100$
and $\omega=2\pi$. In this work, the time unity is years. In order to
interpret these values, we consider $b$ as an annual birth rate, with the
other four rates being annual too. Thus, $\omega=2\pi$ means seasonality with
a period of one year. Regarding the vaccination parameters, the reference
curve (red line) is obtained for $p=v=0$, for the other two $p=0.25$ and
$v=0.2$ in the blue curve, and $v=0.4$ in the green one. We adopt the initial
condition $(s_{0},e_{0},i_{0})=(0.9,0,0.1)$ and discard the first $10^{5}$
integration steps as transient. For these configurations, the solutions of the
system are limit cycles, which are projected onto the plane $i\times s$ in
Fig. 3(a). These cycles contract as the vaccination rate increase, with a
significant decrease in the local maxima of curve $i$, as shown in Fig. 3(b).
We find that the presence of a vaccine campaign causes a reduction of the
maximum number of infects at a time and this effect is amplified with the $v$
increment. Compared to the reference curve, a constant vaccination rate of
one-quarter of the newborn population and one-fifth of the susceptible ones,
for year, leads to approximately a $21.4\%$ reduction in peak infections.
Intensifying vaccination to a rate of $40\%$ of the susceptible population per
year, this reduction reaches $\approx 38.9\%$.

For the next numerical simulations, based on the Gabrick’s work Gabrick2023 ,
we consider the parameters $b=0.02$, $\beta_{0}=270$, $\beta_{1}=0.28$,
$\delta=0.25$, $\alpha=100$, $\gamma=100$ and $\omega=2\pi$. To compute the
impact of the vaccine on the system, over a given finite time interval, we
consider the ratio

$\theta=\frac{A_{\rm v}}{A_{0}},$ (24)

where $A_{\rm v}$ and $A_{0}$ are the areas under the $i$ curve with and
without vaccine, respectively. Being

$\displaystyle A_{\rm v}$ $\displaystyle=\int\limits_{t_{0}}^{t_{0}+\Delta
t}i(t)dt,$ (25) $\displaystyle A_{0}$
$\displaystyle=\int\limits_{t_{0}}^{t_{0}+\Delta t}i(t;p=v=0)dt,$ (26)

where $A_{0}$ is the reference value. From this relation we conclude that if
$\theta<1$, then the vaccination campaign reduces the total number of infected
individuals, and if $\theta>1$, the campaign has the undesired effect of
increasing it. $\theta=1$ means that the vaccination has no effect on the
total infections over time.

Figure 4 shows the influence of vaccination in two different strategies. In
panel (a) we fixed the fraction of newborns $p=0.25$ and get $\theta$ as a
function of the vaccination rate $v\in[0,1]$, with the horizontal axis
discredited in steps of $\Delta v=10^{-2}$. We evolve the system from the same
initial condition used in Fig. 3 and calculate $\theta$ both for the first
$75$ years without transient discard (black curve) and, past $100$ years of
evolution, for the last $75$ years (blue curve). Without discarding the early
years, $\theta$ presents a linear decay as a function of $v$, increasing until
$v\approx 0.44$, and stabilising in $\theta\approx 0.075$. This value is due
to the nonzero area associated with the spread of infection from the initial
condition. In the blue curve we see that the DFE is reached for $v\geq 0.44$
(grey background). Furthermore, it is worth mention that in our simulation if
we consider a large number of years, e.g., 100, 125 years, the DFE also is
reached for $v\geq 0.44$. Figure 4(b) displays $\theta$ as a function of $p$,
with $v=0.1$. This result shows when whole newborns are vaccinated occurs a
reduction $\approx 12.5\%$ of the total infected individuals.

Figure 4:  (a) $\theta$ as function of rate $v$ with $p=0.25$. DFE is reached
for $v=0.44$ (grey background) (b) $\theta$ as function of newborn vaccine $p$
with $v=0.1$. Vaccination of newborns reduced $\approx 12.5\%$ the total of
infected individuals. Horizontal axis discredited in steps of $10^{-2}$.
Results of the first $75$ years without transient discard (black curve) and
evaluated the last $75$ years (blue curve) after $100$ years of the system
evolution. We consider $b=0.02$, $\beta_{0}=270$, $\beta_{1}=0.28$,
$\delta=0.25$, $\alpha=100$, $\gamma=100$ and $\omega=2\pi$; initial condition
$(s_{0},e_{0},i_{0})=(0.9,0,0.1)$.

We evaluate the effects of combining different newborns and susceptible
vaccination rates in the parameter plane $p\times v$, with $\theta$ in colour
scale, displayed in Fig. 5. Our results show the prevalence of $v$ over $p$ in
reducing the spread of infection. This is due to the intrinsic characteristics
of the model, if $\delta$ decreases, the parameter $p$ becomes more relevant.
DFE is achieved within $75$ years from the onset of infection by $(v,p)$ pairs
from the dashed white line to the right of it, defined by the straight line
equation $c_{1}:p+18.55v-8.22=0$. Extinction of infection in the host
population occurs when $e=i=0$ and can be verified in simulations with a given
numerical precision. In order to determine the boundary line $c_{1}$ more
accurately, we integrate the system in the region defined by $v\in[0.39,0.46]$
and $p\in[0,1]$ in a uniform grid of $201\times 201$ points. Still in Fig. 5,
the blue band indicates that the constant vaccination strategy leads to a
significant reduction of the total infections, computed around $20\%$ of the
reference. In the green band, where $v\approx 0.2$, the reduction is around
$50\%$. For small vaccination rates, with $v\approx 10\%$ and less, the amount
of infections decreases only slightly, being around $70\%$ to $80\%$ (gradient
from orange to red) of the result without vaccination, even larger than $90\%$
(magenta band).

Figure 5: Parameter plane $p\times v$ with $\theta$ in colour scale. Both axis
discretised in steps of $10^{-2}$. We consider $b=0.02$, $\beta_{0}=270$,
$\beta_{1}=0.28$, $\delta=0.25$, $\alpha=100$, $\gamma=100$ and $\omega=2\pi$.
$\theta$ is calculated for the first 75 years of infection, without transient
discard and with initial condition $(s_{0},e_{0},i_{0})=(0.9,0,0.1)$. DFE
occurs from the dashed white line ($c_{1}$) to the right of it, where $p\geq
8.22-18.55v$.

The vaccination does not only affects the number of infected individuals, but
also the dynamics. Such effects are examined by hysteresis type bifurcations
diagrams (HTBD), as shown in Fig. 6. This kind of bifurcation diagram is
generated by evolving the system in a given discretized interval of a control
parameter in both directions along the horizontal axis, first in its growth
(red points), then in decrease (blue points), assuming the final state of the
system at the current parameter value as the initial condition for the next
one. This numerical technique is especially useful for finding bi-stability,
when it exists Medeiros2017 . Figure 6(a) displays a HTBD of $i$ local maxima
values ($i_{\rm max}$) as a function of the susceptible vaccination rate.

Figure 6: (a) Hysteresis type bifurcations diagram (HTBD) for $i$ local maxima
values ($i_{\rm max}$) as function of $v\in[0,0.4]$ discretised in steps of
$4\times 10^{-4}$. We consider $b=0.02$, $\beta_{0}=270$, $\beta_{1}=0.28$,
$\delta=0.25$, $\alpha=100$, $\gamma=100$, $\omega=2\pi$ and $p=0.25$. For
each value of $v$, a transient of $10^{5}$ integration steps is discarded and
computed $i_{\rm max}$ over the last $75$ years of the infection. The red
points are in the forward $v$ direction and blue points in the backward
direction. Grey and green backgrounds highlight bi-stability intervals.
Projections in the plane $i\times s$ of periodic (blue line) and chaotic (red
line) attractors coexisting for (b) $v=0.16$ and (c) $v=0.26$.

To compute the peaks in the $i$ time series, we discard the first $10^{5}$
integration steps as transient and evaluate the evolution of the infection
over the last $75$ years in the simulation. For the considered parameters, the
bi-stability dynamics between chaotic and periodic attractors exist without
vaccination Gabrick2023 . The inclusion of vaccination is able to destroy the
bi-stability for some parameters values. However, the bi-stability remains in
two ranges of the control parameter, highlighted in Fig. 6(a), by the grey and
green background columns. In the approximate interval $v\in(0.155,0.165)$
(grey background) the periodic orbit (blue dots) has a maximum value $i_{\rm
max}\approx 0.032$, while the chaotic one (red dots) has a smaller maximum of
$i_{\rm max}\approx 0.026$. Also in the interval $v\in(0.256,0.263)$ (green
background), we observe the periodic orbit with an extreme value of $i_{\rm
max}$ slightly higher than that seen from the chaotic one. Chaotic behaviour
can be related to non-predictability in infectious diseases spread Gabrick2023
. Prior knowledge of these behaviors or how they can be modified is crucial
for epidemic control strategies. Figure 6(b) and 6(c) display the $i\times s$
projection of the attractors obtained for the highlighted values in the bi-
stability for, respectively, $v=0.16$ and $v=0.26$. Both plots show the
periodic solution (blue line) with a higher peak than the chaotic one (red
line). For $v=0.16$, we observe a chaotic attractor from the initial condition
$P_{1}(s_{0},e_{0},i_{0})=(0.9,0,0.1)$, while the periodic one emerges for
$P_{2}(s_{0},e_{0},i_{0})=(0.3,0.03,0.03)$. However, for $v=0.26$, $P_{1}$
leads to the periodic attractor and $P_{2}$ to the chaotic one. The
magnification in the panel (c) exhibits a small amplitude local maximum of
$i$, occurring before the large increase in infectious cases and showing the
two maxima as obtained in the bifurcation diagram.

##  IV Pulsed vaccination strategy

Periodic vaccination consists of the appliance of vaccine in a certain
fraction of the population in periodic time intervals. Examples of this
strategy is considered in seasonal influenza Alkhamis2022  and in control of
childhood viral infections, such as measles and polio Nokes1995 .

Figure 7: Schematic illustration of pulsed $v$ with maximum $v_{\rm max}$
(dashed grey line level), and minimum level $v_{\rm min}$ (solid grey line).
Pulses in four different timings relative to $\beta(t)$ (black curve): at
inflections (green and magenta vertical lines) and at the local minimum and
maximum, blue and red lines, respectively.

In our modelling approach, we consider this immunisation strategy consists of
maintaining a minimum level $v=v_{\rm min}$ during the year and raising the
vaccination rate to $v=v_{\rm max}\geq v_{\rm min}$ in concentrated pulses at
specific times according to the oscillation of the transmissivity. An
illustrative scheme is shown in Fig. 7, where different pulse timings are
superimposed on the periodic curve $\beta(t)$ (black line). We consider four
variants of the same strategy, in which the pulse occurs, relative to the
$\beta(t)$, at every time corresponding to the: inflection point at decreasing
(green line); local minimum (blue line); inflection point at the curve growth
(magenta line); and local maximum (red line). Being the vaccination rate

$\displaystyle v(t;\tau)=\begin{cases}v_{\rm
max},~{}\text{if}~{}~{}\\{t-\tau\\}=0,\\\ v_{\rm
min},~{}\text{otherwise},\end{cases}$ (27)

where $\\{t-\tau\\}$ is the non-integer part of $(t-\tau)$. For the first
variant of the strategy we have $\tau=0.25$, in the second $\tau=0.5$, next
$\tau=0.75$ and in the last one $\tau=0$.

Considering a time interval of $\approx 9$ hours in one day per year, we
investigate the impact on total infected in the first $75$ years of infection,
with $v_{\rm min}=0$ and $v_{\rm max}\in[0,1000]$, where $v=1000$ means that
the entire susceptible population is vaccinated in one pulse, since the
integration step is $10^{-3}$ and the proportion of people vaccinated is given
by the product of the campaign duration with the vaccination rate.
Furthermore, we adopt $p=0$ and the same initial condition used in the
previous section.

Figure 8: $\theta$ as a function of $v_{\rm max}\in[0,1000]$ from four
variants of the pulsed vaccination strategy and obtained in the first $75$
years of infection. Horizontal axis discretised in $1001$ equidistant points.
We consider $b=0.02$, $\beta_{0}=270$, $\beta_{1}=0.28$, $\delta=0.25$,
$\alpha=100$, $\gamma=100$, $\omega=2\pi$ and $p=v_{\rm min}=0$. Initial
condition $(s_{0},e_{0},i_{0})=(0.9,0,0.1)$. Pulse application timing
identified in colours, according to the legend. DFE is achieved (grey
background) within $75$ years for $v_{\rm max}\geq 452$ with
$\tau\in\\{0,0.5\\}$, $v_{\rm max}\geq 459$ for $\tau=0.25$ and $v_{\rm
max}\geq 444$ when $\tau=0.75$.

Figure 8 displays $\theta$ as a function of $v_{\rm max}$, with the horizontal
axis discretised in steps of $\Delta v_{\rm max}=1$, equivalent to increments
of $0.1\%$ in the reach of the immunisation campaign. The four variants of the
strategy are covered, with the pulse applied at the first (green) and second
(magenta) inflections of $\beta(t)$ and, complementary, for both $\tau=0.5$
(blue) and $\tau=0$ (red), at times corresponding of the transmissivity
minimum and maximum, respectively. These results are practically the same.
There is no expressive difference as to when a concentrated immunisation
campaign works, with a slight advantage in applying the pulse at the local
maximum ($\tau=0$) or second inflection ($\tau=0.75$) of the transmissivity
curve. When the campaign coincides with the extremes of the $\beta(t)$ curve
($\tau=0$ and $\tau=0.5$), DFE is achieved within the initial $75$ years (grey
background) from $v_{\rm max}=452$, i.e. vaccinating from $45.2\%$ of the
susceptible population on a single day per year. This threshold is $49.0\%$
for $\tau=0.25$ and $44.4\%$ with $\tau=0.75$, the latter being the best
result.

Figure 9: Parameter plane $v_{\rm max}\times v_{\rm min}$ discretised on a
uniform grid of $101\times 101$ points, with $\theta$ (colour scale)
calculated for the first $75$ years of infection. We consider $b=0.02$,
$\beta_{0}=270$, $\beta_{1}=0.28$, $\delta=0.25$, $\alpha=100$, $\gamma=100$,
$\omega=2\pi$, $p=0$ and $\tau=0.75$. Initial condition
$(s_{0},e_{0},i_{0})=(0.9,0,0.1)$. DFE occurs from the dashed white line
($c_{2}$) to the right of it, where $v_{\rm max}\geq 444-1000v_{\rm min}$.

According to our simulations, considering a not null baseline $v_{\rm min}>0$,
it is possible to reduce the pulse amplitude and yet result in the extinction
of the disease in the host population within $75$ years. In Fig. 9, we show
the parameter plane $v_{\rm max}\times v_{\rm min}$ reinforcing this proposal.
We define $\tau=0.75$, applying the pulse at the moment of inflection in
increasing of the transmissivity, varying $v_{\rm min}\in[0,0.4]$ and $v_{\rm
max}\in[100,400]$ in a uniform grid of $101\times 101$ points. We highlight
the white dashed line $c_{2}:v_{\rm max}+1000v_{\rm min}-444=0$, from which we
obtain DFE. Pairs $(v_{\rm min},v_{\rm min})$ in the blue band bring the total
infected to $\approx 20\%$ of those that would occur without vaccination,
e.g., how we get with the base rate $v_{\rm min}\approx 0.28$ and the
intensive campaign $v_{\rm max}=100$, the latter meaning vaccination of $10\%$
of the susceptible population. Since the base value throughout the year is
around $v_{\rm min}=0.125$, even with $v_{\rm max}=100$, the number of people
infected during the $75$ years simulated is reduced by approximately half
(green band), when compared to the reference case.

Along lines parallel to $c_{2}$, the proportion of susceptible vaccinated
annually is a constant given by

$\rho_{\rm line}=0.001v_{\rm max}+v_{\rm min}.$ (28)

Thus, on line $c_{2}$ we have $\rho_{c_{2}}=0.444$, close to the DFE threshold
value in line $c_{1}$ with $p=0$, shown in Sec. III. In these simulations, an
immunisation campaign that reaches $\approx 22.5\%$ of susceptible annually
reduces to $50\%$ the accumulated of infected people, whereas a rate of
$\approx 38\%$ causes a decrease of $80\%$. Similar values to those obtained
with $v$ constant and without newborns vaccination, displayed in Fig. 5. The
strategy of intensifying immunisation in pulses allows the reduction at
baseline, being $v_{\rm min}$ less than the rates required for the same
results with constant vaccination rate.

##  V Pulsed width vaccination strategy

Extending the proposal of the previous section, we modify the pulsed
vaccination strategy considering, now, a longer duration of the intensified
campaign. Similarly to the formulation of Eq. (27), the pulse is centred at
$t=\tau+k$, $\forall~{}k\in\mathbb{Z}$, but with a width of $D\in[0,1]$, which
corresponds to the duty cycle. Outside the intensive campaign interval,
vaccination assumes a baseline rate $v_{\rm min}$, according to the following
mathematical description

$v(t;\tau)=\begin{cases}v_{\rm
max},~{}\text{if}~{}~{}\left\\{(t-\tau)+\dfrac{D}{2}\right\\}\leq D,\\\ v_{\rm
min},~{}\text{otherwise}.\end{cases}$ (29)

Figure 10: Schematic illustration of pulsed width vaccination. Transmissivity
$\beta(t)$ (black curve) superimposed by vaccination rate $v(t)$ (magenta
curve), it is on two levels: maximum $v_{\rm max}$ (dashed grey line level)
and minimum $v_{\rm min}$. Pulse width $D$ centred on $\tau$, which is a time
instant relative to the seasonal cycle.

Remembering that $\\{x\\}\in[0,1)$ is the non-integer part of the real number
$x$. Figure 10 displays a schematic illustration of this strategy, where the
seasonal transmissivity (black curve) is superimposed on the time-dependent
vaccination rate curve (magenta curve). For instance, in the representation
$D=0.2$, corresponds to one fifth of the year. The concentrated pulse
strategy, according Eq. (27), is recovered when $D=0$. As for $D=1$ we return
to the constant vaccination rate. Based on the best result presented in Sec.
IV, in the next simulations we adopt $\tau=0.75$.

First, we keep $p=0$ and obtain $\theta$ in the parameter planes shown in Fig.
11. As in the previous sections, we take into account the first $75$ years of
infection and consider the same initial condition. Both panels display uniform
grids of $101\times 101$ points.

Figure 11: Parameter planes discretised on a uniform grid of $101\times 101$
points, with $\theta$ (colour scale) calculated for the first $75$ years of
infection. (a) Plane $v_{\rm max}\times D$ with $v_{\rm min}=0$ fixed. DFE
threshold $c_{3}:Dv_{\rm max}=0.44$ and $\theta\approx 0.4$ along the curve
$c_{4}:Dv_{\rm max}=0.28$ (black dashed line). (b) Plane $v_{\rm min}\times D$
with $v_{\rm max}=0.5$ fixed. DFE threshold $c_{5}:v_{\rm
min}=(0.44-0.5D)/(1-D)$ and $\theta\approx 0.4$ along the curve $c_{6}:v_{\rm
min}=(0.28-0.5D)/(1-D)$ (black dashed line). We consider $b=0.02$,
$\beta_{0}=270$, $\beta_{1}=0.28$, $\delta=0.25$, $\alpha=100$, $\gamma=100$,
$\omega=2\pi$, $p=0$ and $\tau=0.75$. Initial condition
$(s_{0},e_{0},i_{0})=(0.9,0,0.1)$. In both panels, DFE occurs from the dashed
white lines ($c_{3}$ and $c_{5}$) to the right, where $\rho\geq 0.44$.

We analyse the results taking into account the annual reaching of the
immunisation campaign, given by

$\rho=v_{\rm min}+D(v_{\rm max}-v_{\rm min}),$ (30)

and the proportion of vaccinated susceptible during each pulse, that is

$\rho_{\rm pulse}=Dv_{\rm max}.$ (31)

In Fig. 11(a) we fix the baseline $v_{\rm min}=0$ and vary the two remaining
parameters of Eq. (29) in the range $0\leq D,v_{\rm max}\leq 1$. In this
configuration $\rho=\rho_{\rm pulse}$ and the DFE threshold is obtained from
$\rho_{\rm pulse}=0.44$, i.e. from the curve $c_{3}:Dv_{\rm max}=0.44$ to the
right of it. Already $\rho_{\rm pulse}=0.28$ leads to $\theta\approx 0.4$
(cyan band), value obtained along the curve $c_{4}:Dv_{\rm max}=0.28$ (dashed
black line). Similar relationships can be obtained for the other $\theta$
values, with the balance between pulse height and duty cycle being the
determining characteristic for reducing the number of infected people over the
$75$ years simulated. Figure 11(b) illustrates the parameter plane $v_{\rm
min}\times D$ with $v_{\rm max}=0.5$, where we evaluate the baseline
vaccination interval $0\leq v_{\rm min}\leq 0.4$. As evidenced in panel (a),
the constant $\theta$ curves are determined by the proportion of susceptible
individuals vaccinated during a year, so we get them from Eq. (30). We check
that DFE is reached from the curve $c_{5}:(0.44-0.5D)/(1-D)$ to the right,
where $44\%$ of susceptible are vaccinated annually. The total infected are
reduced to $\approx 40\%$ around the curve $c_{6}:(0.28-0.5D)/(1-D)$ (black
dashed line), being $28\%$ of susceptible vaccinated along each year.

Inclusion of newborns immunisation changes the minimum proportion of
susceptible immunised annually that leads to DFE. In order to obtain the
boundary surface $c_{\rm DFE}$ as a function of $p$, we use the equation of
the line $c_{1}$, described in Sec. III. Assuming a constant vaccination rate
$v=\rho_{\rm DFE}$ along that line and relating to Eq. (30), we infer

$v_{\rm min}+D(v_{\rm max}-v_{\rm min})-\frac{8.22-p}{18.55}=0.$ (32)

DFE is achieved, during the simulated $75$ years, once

$v_{\rm min}+D(v_{\rm max}-v_{\rm min})\geq\frac{8.22-p}{18.55}.$ (33)

This relation recovers the equations obtained for the thresholds $c_{1}$,
$c_{2}$, $c_{3}$ and $c_{5}$ presented throughout the text. Note that in the
plane, $c_{\rm DFE}$ reduces to curves.

Figure 12: (a) Hysteresis type bifurcations diagram (HTBD) for $i$ local
maxima values ($i_{\rm max}$) as function of $v_{\rm max}\in[0,1]$ discretised
in steps of $10^{-3}$. We consider $b=0.02$, $\beta_{0}=270$,
$\beta_{1}=0.28$, $\delta=0.25$, $\alpha=100$, $\gamma=100$, $\omega=2\pi$,
$p=v_{\rm min}=0$ and $D=0.4$. For each value of $v_{\rm max}$, a transient of
$10^{5}$ integration steps is discarded and computed $i_{\rm max}$ over the
last $75$ years of the infection. Red points are in the forward $v_{\rm max}$
direction and blue ones in the backward direction. Grey and green backgrounds
highlight bi-stability intervals. Projections in the plane $i\times s$ of
periodic (blue line) and chaotic (red line) attractors coexisting for (b)
$v=0.01$ and (c) $v=0.56$.

In addition, the bi-stability depends on the pulsed vaccination parameters.
Figure 12(a) displays a HTBD of infected local maxima values ($i_{\rm max}$)
as a function of $v_{\rm max}\in[0,1]$, with $p=v_{\rm min}=0$ and $D=0.4$.
The results yielded in forward direction are red points and backward ones are
the blue points. The peaks are computed discarding the first $10^{5}$
integration steps as transient and evaluate the $i$ series over the last $75$
years in the simulation. Similarly to the case of a constant vaccination rate,
shown in Fig. 6(a), there are two ranges (highlighted backgrounds) of the
control parameter where periodic and chaotic orbits coexist. For $v_{\rm
max}\in[0,0.014)$ (grey background) the periodic orbit (blue dots) has a
maximum value $i_{\rm max}\approx 0.049$, However, the chaotic one (red dots)
has a smaller maximum of $i_{\rm max}\approx 0.036$, as also showed by the
attractors projection in Fig. 12(b) for $v_{\rm max}=0.01$. In the next
interval $v_{\rm max}\in(0.542,0.577)$ (green background), the periodic orbit
presents a maximum value $i_{\rm max}\approx 0.019$, while the chaotic one has
an extreme $\approx 0.015$, i.e., the periodic solution has peak infections
values $26.7\%$ above the chaotic case. This difference between the peak of
infectious is better observed in the $i\times s$ projection of the attractors,
as displayed in Fig. 12(c), for $v_{\rm max}=0.56$. The chaotic solutions (red
line) in both Fig. 12(b) and 12(c) are obtained from the initial condition
$P_{1}$, the periodic attractors (blue line) are from (b)
$P_{3}(s_{0},e_{0},i_{0})=(0.27,0.05,0.05)$ and (c)
$P_{4}(s_{0},e_{0},i_{0})=(0.28,0.02,0.02)$.

##  VI Conclusions

In this work, we study a SEIRS model with seasonal transmissivity and
vaccination control. Considering an autonomous version of the system, we
obtain two equilibrium points: the DFE and the endemic fixed point. For the
DFE, we established a limit for the vaccination rate as function of the other
parameters in the model. We find that the endemic solution can only exist if
the DFE is not stable. By numerical simulations, we explore three different
vaccination strategies in the non-autonomous case: constant, pulsed and pulsed
width vaccination strategy. All are based on immunising a proportion of
susceptible individuals, as well as a fraction of newborns. In order to
evaluate these three different strategies, we analyse the accumulated infected
over a simulated interval of $75$ years:

$i$) The first one refers to a constant vaccination rate applied to the
susceptible individuals. Our results show that the immunisation of newborns is
able to slightly reduce the total number of infected in the simulated time
interval. However, vaccination of susceptible is more efficient, with a
constant rate of $\approx 44\%$, or greater, leading to the extinction of the
infection in the host population. We also identify the occurrence of bi-
stability as a function of the vaccination rate, with periodic orbits
presenting values of infected higher than the chaotic ones.

$ii$) The second strategy is the pulsed vaccination, in which the immunisation
campaign operates intensively on a single day per year, more precisely and
according to the simulations, just $9$ hours on one day in each year. The
baseline vaccination rate is a small constant. Pulses are applied annually at
timings given in relation to seasonality, thus modelling the campaign acting
according to the transmissivity variation. The timing of the immunisation
pulse does not lead to a significant difference in the accumulated number of
infected. Even so, acting at the inflection point of the rise in the
transmissivity curve proved to be slightly more effective, reducing the
percentage of susceptibles vaccinated during the pulse to result in DFE. The
baseline and intense immunisation rate has a linear relationship with the
reduction of disease spread. Hence, the total number of cases depends on the
proportion of susceptible vaccinated annually, and not on independent
immunisation rates.

$iii$) In our last strategy, we applied a pulsed width vaccination campaign.
We discover that the disease spread depends on a non-linear relationship
between vaccination rates and the duty cycle of the campaign. We derive the
relationship between the immunisation parameters that leads to the extinction
of the infection within the simulated $75$ years. It depends on the baseline
and intensive vaccination rates, the campaign duty cycle and the proportion of
immunised newborns. We find that the reduction of the accumulated infected
number depends on the annual vaccination rate of the susceptible population,
this finding being valid for all three strategies, mainly in the first one.
For example, regardless of the strategy, given a vaccination of $\approx 38\%$
of susceptible each year, the cumulative case of infection is reduced to only
$20\%$ of what it would be in a no-vaccination situation. Since we study the
model with normalised variables, we are always dealing with proportions of the
total population. For the constant population approximation, a direct
transition to the number of infected people over time is valid, whereas for
the more general case, where the population may vary, we must take the results
into account as proportions.

We plan to study, in future works, the effects of periodic and perturbed
vaccination in chaotic dynamics and bi-stable parameter rages. In addition, we
plan to fill some points opened in this work, which is exploring the effects
of stochastic process in vaccination campaign.

## Acknowledgements

The authors thank the financial support from the Brazilian Federal Agencies
(CNPq); São Paulo Research Foundation (FAPESP) under Grant Nos. 2021/12232-0,
2018/03211-6, 2022/13761-9; Coordenação de Aperfeiçoamento de Pessoal de Nível
Superior (CAPES); Fundação Araucária. R.L.V. received partial financial
support from the following Brazilian government agencies: CNPq (403120/2021-7,
301019/2019-3), CAPES (88881.143103/2017-01), FAPESP (2022/04251-7). E.C.G.
received partial financial support from Coordenação de Aperfeiçoamento de
Pessoal de Nível Superior - Brasil (CAPES) - Finance Code
88881.846051/2023-01. We thank 105 Group Science (www.105groupscience.com).

## References

## References

  * (1) Xia, S., Gullickson, C.C., Metcalf, C.J.E., Grenfell, B.T., Mina, M.J., 2023. Assessing the Effects of Measles Virus Infections on Childhood Infectious Disease Mortality in Brazil. The Journal of Infectious Diseases 227,133-40.
  * (2) Glatter, K.A., Finkelman, P., 2021. History of the Plague: An ancient pandemic for the Age of COVID-19. The American Journal of Medicine 134(2), 176-181.
  * (3) Tumpey, T.M., Basler, C.F., Aguilar, P.V., Zeng, H., Solórzano, A., Swayne, D.E., Cox, N.J, Katz, J.M., Taubenberger, J.K., Palese, P., García-Sastre, A., 2005. Spanish influenza pandemic virus. Science 310, 77–80.
  * (4) Manchein, C., Brugnago, E.L., da Silva, R.M., Mendes, C.F.O., Beims, M.W., 2020\. Strong correlations between power-law growth of COVID-19 in four continents and the inefficiency of soft quarantine strategies. Chaos 30, 041102.
  * (5) Aguiar, M., Ballesteros, S., Kooi, B.W., Stollenwerk, N., 2011. The role of seasonality and import in a minimalistic multi-strain dengue model capturing differences between primary and secondary infections: Complex dynamics and its implications for data analysis. Journal of Theoretical Biology 289, 181-196.
  * (6) Dalal, N., Greenhalgh, D., Mao, X., 2008. A stochastic model for internal HIV dynamics. Journal of Mathematical Analysis and Applications 341, 1084-1101.
  * (7) Meyers, L.A., Pourbohloul, B., Newman, M.E.J., Skowronski, D.M., Brunham,R.C., 2005\. Network theory and SARS: predicting outbreak diversity. Journal of Theoretical Biology 232(1), 71-81.
  * (8) Mao, L., Bian, L., 2010. Spatial–temporal transmission of influenza and its health risks in an urbanized area. Computers, Environment and Urban Systems 34(3), 204-215.
  * (9) Scarpino, S.V., Petri, G., 2019. On the predictability of infectious disease outbreaks. Nature Communications 10, 898.
  * (10) Altizer, S., Dobson, A., Hosseini, P., Hudson, P., Pascual, M., Rohani, P., 2006\. Seasonality and the dynamics of infectious diseases. Ecology Letters 9, 467-484.
  * (11) Buonomo, B., Chitnis, N., d’Onofrio, A., 2018. Seasonality in epidemic models: a literature review. Ricerche di Matematica 67, 7-25.
  * (12) Hoshen, M.B., Morse, A.P., 2004. A weather-driven model of malaria transmission. Malaria Journal 3(32), 10.1186/1475-2875-3-32.
  * (13) Finkenstadt, B.F., Grenfell, B.T., 2000. Time series modelling of childhood diseases: a dynamical systems approach. Journal of the Royal Statistical Society Series C Applied Statistics 49(2), 187-205.
  * (14) Grassly, N.C., Fraser, C., 2006. Seasonal infectious disease epidemiology. Proceedings of the Royal Society B: Biological Sciences 273(1600), 2541–2550.
  * (15) Azimaqin, N., Pen, Z., Ren, X., Wei, Y., Liu, X., 2022. Vaccine failure, seasonality and demographic changes associate with mumps outbreaks in Jiangsu Province, China: Age-structured mathematical modelling study. Journal of Theoretical Biology 544, 111125.
  * (16) Aguiar, M., Stollenwerk, N., Kooi, B.W., 2009. Torus bifurcations, isolas and chaotic attractors in a simple dengue fever model with ADE and temporary cross immunity. International Journal of Computer Mathematics 86(10-11), 1867-1877.
  * (17) Mugnaine, M., Gabrick, E.C., Protachevicz, P.R., Iarosz, K.C., de Souza, S.L.T., Almeida, A.C.L., Batista, A.M., Caldas, I.L., Szezech Jr., J.D., Viana, R.L., 2022. Control attenuation and temporary immunity in a cellular automata SEIR epidemic model. Chaos, Solitons and Fractals 155, 111784.
  * (18) Batista, A.M., de Souza, S.L.T., Iarosz, K.C., Almeida, A.C.L., Szezech Jr., J.D., Gabrick, E.C., Mugnaine, M., dos Santos, G.L., Caldas, I.L., 2021. Simulation of deterministic compartmental models for infectious diseases dynamics. Revista Brasileira de Ensino de Física 43, e20210171.
  * (19) Grenfell, B.T., Bolker, B.M., Kleczkowski, A., 1995. Seasonality and Extinction in Chaotic Metapopulations. Proceedings: Biological Sciences 259(1354), 97-103.
  * (20) Keeling, M.J., Rohani, P., 2008. Modeling Infectious Diseases in Humans and Animals. Princeton University Press.
  * (21) Bilal, S., Singh, B.K., Prasad, A., Michael, E., 2016. Effects of quasiperiodic forcing in epidemic models. Chaos 26, 093115.
  * (22) Tanaka, G., Aihara, K., 2013. Effects of seasonal variation patterns on recurrent outbreaks in epidemic models. Journal of Theoretical Biology 317, 87-95.
  * (23) Olsen, L.F., Schaffer, W.M., 1990. Chaos versus noisy periodicity: alternative hypotheses for childhood epidemics. Science 249, 499-504.
  * (24) Gabrick, E.C., Sayari, E., Protachevicz, P.R., Szezech Jr., J.D., Iarosz, K.C., de Souza, S.L.T., Almeida, A.C.L., Viana, R.L., Caldas, I.L., Batista, A.M., 2023. Unpredictability in seasonal infectious diseases spread. Chaos, Solitons and Fractals 166, 113001.
  * (25) London, W.P, Yorke, J.A., 1973. Recurrent outbreaks of measles, chickenpox and mumps: I. seasonal variation in contact rates. American Journal of Epidemiology 98(6), 453–468.
  * (26) Greenhalgh, D., and Moneim, I.A., 2003. SIRS epidemic model and simulations using different types of seasonal contact rate. Systems Analysis Modelling Simulation 43(5), 573-600.
  * (27) de Carvalho, J.P.M., Rodrigues, A.A., 2022. Strange attractors in a dynamical system inspired by a seasonally forced SIR model. Physica D 434, 133268\.
  * (28) Metcalf, C.J.E., Bjornstad, O.N., Grenfell, B.T., Andreasen, V., 2009. Seasonality and comparative dynamics of six childhood infections in pre-vaccination Copenhagen. Proceedings of The Royal Society B 276, 4111–4118.
  * (29) Keeling, M.J., Rohani, P., Grenfell, B.T., 2001. Seasonally forced disease dynamics explored as switching between attractors. Physica D 148, 317–335.
  * (30) Cooper, I., Mondal, A., Antonopoulos, C.G., 2020. A SIR model assumption for the spread of COVID-19 in different communities. Chaos, Solitons and Fractals 139, 110057.
  * (31) Brugnago, E.L., da Silva, R.M., Manchein, C., Beims, M.W., 2020, How relevant is the decision of containment measures against COVID-19 applied ahead of time? Chaos, Solitons and Fractals 140, 110164.
  * (32) de Souza, S.L.T., Batista, A.M., Caldas, I.L., Iarosz, K.C., Szezech Jr., J.D., 2021\. Dynamics of epidemics: Impact of easing restrictions and control of infection spread. Chaos, Solitons and Fractals 142, 110431.
  * (33) Ansari, S., Anvari, M., Pfeffer, O., Molkenthin, N., Moosavi, M.R., Hellmann, F., Heitzig, J., Kurths, J., 2021. Moving the epidemic tipping point through topologically targeted social distancing. The European Physical Journal Special Topics 230, 3273–3280.
  * (34) Mello, B.A., 2020. One-Way Pedestrian Traffic Is a Means of Reducing Personal Encounters in Epidemics. Frontiers in Physics 8, 376.
  * (35) Balsa, C., Lopes, I., Guarda, T., Rufino, J., 2021. Computational simulation of the COVID-19 epidemic with the SEIR stochastic model. Computational and Mathematical Organization Theory 2021, 1-19.
  * (36) Zou, Y., Yang, W., Lai, J., How, J., Lin, W., 2022. Vaccination and quarantine effect on COVID-19 transmission dynamics incorporating Chinese-Spring-Festival travel rush: modeling and simulations. Bulletin of Mathematical Biology 84(30), 1-19.
  * (37) Etxeberria-Etxaniz, M., Alonso-Quesada, S., De la Sen, M., 2020. On an SEIR Epidemic Model with Vaccination of Newborns and Periodic Impulsive Vaccination with Eventual On-Line Adapted Vaccination Strategies to the Varying Levels of the Susceptible Subpopulation. Applied Sciences 10(22), 8296\.
  * (38) De la Sen, M., Alonso-Quesada, S., Ibeas, A., Nistal, R., 2021. On a Discrete SEIR Epidemic Model with Two-Doses Delayed Feedback Vaccination Control on the Susceptible. Vaccines 9(4), 398.
  * (39) Voyseym, M., Clemens, S.A., Madhi, S., Weckx, L., Folegatti, P., Aley, P., Angus, B., Baillie, V., Barnabas, S., Bhorat, Q., Bibi, S., Briner, C., Cicconi, P., Clutterbuck, E., Collins, A., Cutland, C., Darton, T., Dheda, K., Chritina, C., 2021. Single-dose administration and the influence of the timing of the booster dose on immunogenicity and efficacy of ChAdOx1 nCoV-19 (AZD1222) vaccine: a pooled analysis of four randomised trials. The Lancet 397, 881-891.
  * (40) Thompson, H.A., Hogan, A.B., Walker, P.G.T., Winskill, P., Zongo, I., Sagara, I., Tinto, H., Ouedraogo, J.B., Dicko, A., Chandramohan, D., Greenwood, B., Cairns, M., Ghani, A.C., 2022. Seasonal use case for the RTS,S/AS01 malaria vaccine: a mathematical modelling study. Lancet Global Health 10, e1782–92.
  * (41) Molinari, N.A.M., Ortega-Sanchez, I.R., Messonnier, M.L., Thompson, W.W., Wortley, P.M., Weintraub, E., Bridges, C.B., 2007. The annual impact of seasonal influenza in the US: Measuring disease burden and costs. Vaccine 25(27), 5086-5096.
  * (42) Keeling, M.J., Grenfell, B.T., 2002. Understanding the persistence of measles: reconciling theory, simulation and observation. Proceedings of the Royal Society of London. Series B: Biological Sciences 269(1489), 335–343.
  * (43) Gao, S., Liu, Y., Nieto, J.J., Andrade, H., 2011. Seasonality and mixed vaccination strategy in an epidemic model with vertical transmission. Mathematics and Computers in Simulation 81, 1855–1868.
  * (44) Atchison, C., Lopman, B., Edmunds, W.J., 2010. Modelling the seasonality of rotavirus disease and the impact of vaccination in England and Wales. Vaccine 28, 3118–3126.
  * (45) Shulgin, B., Stone, L., Agur, Z., 1998. Pulse vaccination strategy in the SIR epidemic model. Bulletin of Mathematical Biology 60, 1123-1148.
  * (46) Gabrick, E.C., Protachevicz, P.R., Batista, A.M., Iarosz, K.C., de Souza, S.L.T., Almeida, A.C.L., Szezech Jr., J.D., Mugnaine, M., Caldas, I.L., 2022. Effect of two vaccine doses in the SEIR epidemic model using a stochastic cellular automaton. Physica A 597, 127258.
  * (47) Duarte, J., Januário, C., Martins, N., Seoane, J.M., Sanjuán, M.A.F., 2021\. Controlling Infectious Diseases: The Decisive Phase Effect on a Seasonal Vaccination Strategy. International Journal of Bifurcation and Chaos 31(15), 2130044.
  * (48) Wang, X., Pen, H., Shi, B., Jiang, D., Zhang, S., Chen, B., 2019. Optimal vaccination strategy of a constrained time-varying SEIR epidemic model. Communications in nonlinear science and numerical simulations 67, 37-48.
  * (49) Bai, Z., Zhou, Y., 2012. Global dynamics of an SEIRS epidemic model with periodic vaccination and seasonal contact rate. Nonlinear Analysis: Real World Applications 13, 1060–1068.
  * (50) Moneim, I.A., Greenhalgh, D., 2005. Use of a periodic vaccination strategy to control the spread of epidemics with seasonally varying contact rate. Mathematical Biosciences and Engineering 2(3), 591–611.
  * (51) Grabenstein, J.D., Nevin, R.L. 2006. Mass Immunization Programs: Principles and Standards. CTMI 304, 31-51.
  * (52) Alkhamis, A.K., Hosny, M. 2022. A Synthesis of Pulse Influenza Vaccination Policies Using an Efficient Controlled Elitism Non-Dominated Sorting Genetic Algorithm (CENSGA). Electronics 11(22), 3711.
  * (53) Song, P., Lou, Y., Xiao, Y., 2019. A spatial SEIRS reaction–diffusion model in heterogeneous environment. Journal of Differential Equations 267, 5084-5114.
  * (54) Cooke, K.L., van den Driessche, P., 1996. Analysis of an SEIRS epidemic model with two delays. Journal of Mathematical Biology 35, 240-260.
  * (55) Ho, S.H., He, D., Eftimie, R., 2019. Mathematical models of transmission dynamics and vaccine strategies in Hong Kong during the 2017–2018 winter influenza season. Journal of Theoretical Biology 476, 74-94.
  * (56) Metcalf, C.J.E., Lessler, J., Klepac, P., Cutts, F., Grenfell, B.T., 2012. Impact of birth rate, seasonality and transmission rate on minimum levels of coverage needed for rubella vaccination. Epidemiology and Infection 140(12), 2290-2301.
  * (57) Nokes, D.J., Swinton, J., 1995. The control of childhood viral infections by pulse vaccination. Journal of Mathematics Applied in Medicine and Biology 12, 29-53.
  * (58) Agur, Z., Cojocaru, L., Mazor, G., Anderson, R.M., Danon, Y.L., 1993. Pulse mass measles vaccination across age cohorts. Proceedings of the National Academy of Sciences of the USA 90, 11698-11702.
  * (59) Bjørnstad, O.N., 2018. Epidemics: Models and Data using R. Springer.
  * (60) Hernández Guillén, J.D., Martín del Rey, A., Hernández Encinas, L., 2017. Study of the stability of a SEIRS model for computer worm propagation. Physica A: Statistical Mechanics and its Applications 479, 411-421.
  * (61) Tél, T., Gruiz, M., 2006. Chaotic Dynamics: An introduction based on Classical Mechanics. Cambridge University Press.
  * (62) Alligood, K. T., Sauer, T. D., Yorke, J. A., 1996. Chaos: an introduction to dynamical systems. Springer.
  * (63) Mclean, A., Anderson, R. 1988. Measles in developing countries. Part II. The predicted impact of mass vaccination. Epidemiology and Infection 100(3), 419-442.
  * (64) Fenner, F. Henderson, D.A., Arita, I., Jezek, Z., Ladnyi, I.D., World Health Organization. 1988. Smallpox and its eradication. History of international public health 6, 1371-1409.
  * (65) Bagcchi, S. 2021. The world’s largest COVID-19 vaccination campaign. The Lancet Infectious Diseases 21(3), 323.
  * (66) Anderson, R.M., May, R.M., 1992. Infectious diseases of humans: dynamics and control. Oxford University Press.
  * (67) Wesley, C.L., Allen, L.J.S. 2009. The basic reproduction number in epidemic models with periodic demographics. Journal of Biological Dynamics 3(2-3), 116-129.
  * (68) Meinsma, G. 1995, Elementary proof of the Routh-Hurwitz test. Systems & Control Letters 25(4), 237-242.
  * (69) Butcher, J.C., 1987. The Numerical Analysis of Ordinary Differential Equations: Runge-Kutta and General Linear Methods. Wiley, Wiley-Interscience.
  * (70) Medeiros, E.S., Caldas, I.L., Baptista, M.S., Feudel, U., 2017. Trapping phenomenon attenuates the consequences of Tipping Points for limit cycles. Scientific Reports 7, 42351.