\documentclass[a4paper,10.80pt,twoside]{article}
\usepackage{informat}
\usepackage{epsfig}

\usepackage{amsmath}
\usepackage{amsfonts} % for Natural numbers
\usepackage[svgnames]{xcolor}
\usepackage{listings}
\usepackage{comment} 
\usepackage{subfig}
\usepackage{amsthm}
\usepackage[utf8]{inputenc}
\usepackage[greek,english]{babel}
\usepackage{tabularx}
\usepackage{array}

\lstset{language=R,
	frame=tb,
	basicstyle=\small\ttfamily,
	stringstyle=\color{DarkGreen},
	otherkeywords={0,1,2,3,4,5,6,7,8,9},
	morekeywords={TRUE,FALSE},
	deletekeywords={data,frame,length,as,character},
	keywordstyle=\color{blue},
	numberstyle=\tiny\color{gray},
	tabsize=4
}

\theoremstyle{definition}
\newtheorem{definition}{Definition}[section]
\newtheorem{theorem}{Theorem}
\usepackage{float}
\usepackage{blindtext}
\usepackage{hyperref}
\usepackage{multicol} %ORIGINAL
\usepackage{parcolumns}
%\usepackage{physics}
%\usepackage{array,amsmath,booktabs}


%\addbibresource{Paper_format.bib}
%\usepackage[backend=bibtex]{biblatex}

%\usepackage{fontspec}
%\setmainfont{Times New Roman}

\begin{document}
%\setcounter{page}{555}
  \title{Neural Networks and Sigmoidal Models for the Waves of the Coronavirus}
   \author{Oliver Amadeo Vilca Huayta \\
      Departamento de Ingenier�a de Sistemas, Universidad Nacional del Altiplano, Peru\\
      ovilca@unap.edu.pe}
  \titleodd{Neural Networks and Sigmoidal Models for the Waves of the Coronavirus}
  \authoreven{O. Vilca}
  \keywords{Brazil, COVID-19, death number, Italy, logistic sigmoidal model, pandemic, Peru, Switzerland, artificial neural networks}
  \received{November 1, 2024}

  \abstract{The rapid spread of the coronavirus disease 2019 (COVID-19) pandemic has caused enormous problems---mainly a large number of deaths---and, so, it is important to construct epidemiological models for forecasting and prevention.
  	The main objective of this study is to develop a novel and efficient model $F(x)$, treating the cumulative number of deaths due to COVID-19 using two competing approaches: the logistic sigmoidal function and an artificial neural network. We also aim to estimate models for the death rate.
  	This research is retrospective and longitudinal.
  	Data were downloaded from Johns Hopkins University for four countries: Italy, Mexico, Brazil, Switzerland, and Peru.
  	It is shown that the basic logistic function is optimal---that is, it has a running time that cannot be improved---and a new variant of the function is obtained and used. 
  	The novel contribution of this work is the general model that can be used on data that are not necessarily epidemiological; and with any function, where, as $x$ approaches positive or negative infinity, the function tends to a constant value. Also, the method of its construction, and the calculation of high-order derivatives, which allow for the development of an elegant and useful model, as well as justification of the term ``Wave(s)''. The inflection points (i.e., the place where the concavity changes) are also obtained and validated. Finally, they were applied to real-world datasets. 
  	It is concluded that the obtained models have a high fit, with acceptable root mean squared error (RMSE) and Pearson correlation coefficients greater than 0.999. The models are representative, predictive, and optimal and, so, useful to achieve a better understanding of the pandemic and improve future public health response.}

\abstractSi{}

\maketitle

\section{Introduction}  \label{sec:Intro}
%Context and significance
The coronavirus disease 2019 (COVID-19) is highly contagious disease~\cite{sajed_corona_2020,din_mathematical_2020} that has caused a high number of deaths and continues to damage the economy and society, to name just a few of its impacts. 
%The World Health Organization (WHO) has been reporting associated deaths and infections for the year 2024~\cite{WHOcovid-19}.
On October 4, 2024, the Heads of the International Monetary Fund (IMF), the World Bank Group (WBG), and the World Health Organization (WHO) have agreed on broad principles for cooperation on pandemic preparedness (\href{https://www.imf.org/en/News/Articles/2024/10/04/pr24346-imf-wbg-who-step-up-cooperation-on-pandemic-preparedness}{https://www.imf.org}; accessed on 7 October 2024).
The development and study of epidemiological models can be essential to prevent, predict, or mitigate future epidemics or pandemics~\cite{ogunleye_response_2020,bhagat_novel_2022,tabish_covid-19_2020}.
They are also necessary for decision making; for example, allowing public health institutions to take reliable measures in order to preserve people's lives~\cite{kumar_monitoring_2020,alsalem_multi-criteria_2022}. In addition, comparison between countries is considered essential for the control of COVID-19~\cite{pearce_comparisons_2020}.

In this study, four countries with different characteristics were selected for the purposes of modeling, comparison, and research, namely,
Italy, Mexico, Brazil, Switzerland, and Peru, for which COVID-19 has resulted in a total of 188,322, 699,276, 333,188, 14,210, and 219,539 deaths, respectively, according to the accumulated death time-series of Johns Hopkins University (https://coronavirus.jhu.edu; accessed on 1 January
2024), without considering the deaths registered after the end of the pandemic. The number of deaths and mortality rates vary between countries~\cite{jawad_hashim_population_2020} and, in this work, we demonstrate that they follow a strict sigmoidal behavior typical of epidemics or pandemics.

%Method
The Gompertz model, due to the definition of its formula $f(x) = a\cdot e^{-e^{b-c\cdot x}}$ (i.e., double exponential), requires more elementary operations than the most basic of the sigmoidal functions. Therefore, in this research, we demonstrate that the logistic function is the most efficient and optimal for numerical algorithms (e.g., machine learning and non-linear regression). Therefore, given the structure and functionality of the logistic function, we recommend its use in future work.

Additionally, we adapted and obtained a new variant of the logistic model, which we call logistic 2 or just the logistic model $g(x)=\frac{H}{1+\mathrm{2}^{A-Bx}}$, which does not lose any of the original model's properties or qualities. Furthermore, both functions are equivalent.

One of the major contributions of this work is the general model and the method of its construction. Each function $g(x)_i$ corresponds to wave $i$. It has the following characteristics or advantages: It does not require any form of integration or union of functions, does not need preliminary work, and it does not require to add a dummy variable to the data (or to the function) as in previous works~\cite{Vilca-Huayta2022,informatics11020040}. Mainly for this reason and the previous ones it is more efficient:

\begin{eqnarray} \label{sec:FFinal1}
	F \left( x \right) &=& \frac{H_1}{1+\mathrm{2}^{A_1 - B_1 x}} + \frac{H_2}{1+\mathrm{2}^{A_2 - B_2 x}} + \cdots + \nonumber\\ 
	&& \frac{H_i}{1+\mathrm{2}^{A_i - B_i x}} + \cdots + \frac{H_n}{1+\mathrm{2}^{A_n - B_n x}} \nonumber
\end{eqnarray}

where $F(x)$ is the expected cumulative number of deaths, $x$ is the number of days since the first case, $n$ is the number of waves, $1 \leq i\leq n$, $H_i$ is the height of the wave $i$, $A_i$ and $B_i$ are constants related to the pandemic.


%Main aim of the work
The main objective of this study is to develop a novel and efficient model $F(x)$, treating the cumulative number of deaths due to COVID-19 since the pandemic began as a variable, using the logistic function. Furthermore, we compare the proposed model with an artificial neural network (ANN) which, in this case, constitutes an ideal opponent.

%Results
The first derivative is calculated to explain the death rate and the magnitude of the COVID-19 waves, where even its graph represents the considered characteristics. In this way, justification is provided for the denomination of COVID-19 ``waves.'' Next, the second derivative is used to calculate the vertical inflection points (i.e., the point where the sigmoidal function changes from concave upward to concave downward). At this point, only one last issue remains to be resolved: the second derivative can report false inflection points. Therefore, the inflection points are validated through the use of the third derivative. The use of high-order derivatives---specifically, the first three derivatives---allows for an elegant and irrefutable mathematical analysis of any model in general. Then, the proposed model is applied to four representative countries, yielding novel results.

In ANN and the proposed model, the differences were minimal. The disadvantage of the ANN is a considerably longer calculation time, when compared to sigmoidal functions.

This paper is structured as follows: Section \ref{sec:Related} describes related works, Section \ref{sec:MaMethods} presents the materials and methods, and describes the used data. Section \ref{sec:Results} demonstrates the results, including those related to the ANN model, the coupled model, high-order derivatives, and the case study. Finally, the conclusions and directions for future work are presented in Section \ref{sec:conclusionAndFutureWork}.

\begin{table*}[h]
	\begin{center}
		\begin{tabularx}{\textwidth}{ >{\raggedright\arraybackslash}X | >{\centering\arraybackslash}X }
			\hline
			\textbf{In this work} & \textbf{In Key publication} \\
			\hline
			It is not necessary to estimate any horizontal inflection point as a preliminary procedure to obtain $F(x)$. & It is necessary to calculate the horizontal inflection points.  \\
			\hline
			The logistic sigmoidal function was used. It does not require the dummy variable $q$.        & The Boltzmann sigmoidal function was used. It requires the insertion of a dummy variable $q$ in the model and in the database.\\
			\hline
			A general model $F$ was obtained, whose elements are extended to functions that are not necessarily sigmoidal, that is, with $\lim_{x \to - \infty}h(x) = L$ and $\lim_{x \to \infty}h(x) = R$, where $L$ and $R$ are constants. Finally, it was demonstrated with mathematical induction. It was applied to four case studies, where the Peruvian model is:
			\newline
			$ F(x) = \frac{85.94}{1+2^{6.07-0.048*x }}  +  \frac{114.1}{1+2^{15.96-0.041*x }}  +  \frac{13.22}{1+2^{44.35-0.063*x }}  + 
			\frac{3.4}{1+2^{73.16-0.082*x }}  +  \frac{3.33}{1+2^{46.8-0.045*x }} $    & The Peruvian model is: \newline $F(x,q) = \frac{(2-q)(3-q)(4-q)(5-q)}{24}\left\{\frac{88.05}{1+e^\frac{\left(128.95-x\right)}{30.88}}\right\}-\frac{(1-q)(3-q)(4-q)(5-q)}{6}\left\{\frac{113.13}{1+e^\frac{\left(390.65-x\right)}{34.36}}+86.49\right\}+\frac{(1-q)(2-q)(4-q)(5-q)}{4}\left\{\frac{14.18}{1+e^\frac{\left(698.06-x\right)}{26.51}}+199.34\right\}-\frac{(1-q)(2-q)(3-q)(5-q)}{6}\left\{\frac{3.61}{1+e^\frac{\left(898.54-x\right)}{17.99}}+213.37\right\}+\frac{(1-q)(2-q)(3-q)(4-q)}{24}\left\{\frac{2.84}{1+e^\frac{\left(1039.82-x\right)}{25.42}}+216.95\right\} $ \\
			\hline			
			$\sigma = 0.999905$ (the largest). The vertical inflection points and higher order derivatives were calculated. & $\sigma =  0.999$. The vertical inflection points and higher order derivatives were not calculated.\\			
			%\hline
			%\multicolumn{2}{c}{In relation to artificial neural networks}  \\
			%\hline
			%\hline
			%The R programming language was used. The RMSE was used and additionally the Pearson correlation coefficient. & The Python programming language was used. The Pearson correlation coefficient was used. \\
			%\hline
			%The models were applied in four case studies. And a comparative study was carried out. & The model was applied in only one case study (Peru).\\
			\hline
		\end{tabularx}
	\end{center}
	\caption{Main differences in relation to the key publication}
	\label{table:Diffes}
\end{table*}

\section{Related works} \label{sec:Related}

%The current state of the research
In relation to key publications: While one study has focused on the Peruvian case, it was limited to the calculation of a function with the Gompertz model and did not apply ANNs, and was carried out with incomplete data as it was conducted during the pandemic period~\cite{Vilca-Huayta2022}. 
There is also a study in which $F(x)$ was calculated using the Boltzmann sigmoidal function for a case study (Peru), which was correlated with the social isolation measures in Peru (qualitative variable)~\cite{informatics11020040}, Table~\ref{table:Diffes} illustrates the main differences.


Regarding ANN, several works have used artificial neural networks for prediction. In India,  research has been carried out on the Multilayer perceptron~\cite{sujathmachine2020}. In Brazil, a long short-term memory (LSTM) model has been used~\cite{fernandeslong2022}, which obtained a correlation coefficient of 0.96. In the United States and India, various artificial neural networks have been used with data on confirmed infections and deaths from COVID-19, of which the convolutional LSTM stood out, which had greater precision and a lower error~\cite{SHASTRI2020110227}. Additionally, in China and other countries, it has been concluded that ANNs are adequate to predict global infections and deaths due to COVID-19~\cite{namasudra_nonlinear_2023}.

Regarding the use of the logistic sigmoidal function to analyze pandemics, there have not been many related studies. One study used the extended logistic function in this context~\cite{abolmaali_comparative_2021}; however, the authors did not compare it with an ANN, and did not calculate derivatives or inflection points. 


\section{Materials and Methods} \label{sec:MaMethods}

It should be noted that this study consider other model, another data, and use high-order derivatives of $F$. And new simple way to couple functions $g(x)$.

\subsection{Data Set}
The database of the total number of deaths used in the current study is publicly available (open access) from the Coronavirus Disease Data Repository of Johns Hopkins University 
(\href{https://github.com/CSSEGISandData/COVID-19}{https://github.com/CSSEGISandData/COVID-19}; accessed on 1 January 2024). Complete data for different countries are also available~\cite{DatasetCovid19}.

In order to achieve good results, it is necessary to work with complete data from the beginning to the end of the pandemic. Studying partial data may not be sufficient. In addition, modeling a single wave~\cite{haouari_particle_2021} can be relatively easy. In this work, no function integration method is required~\cite{Vilca-Huayta2022}, and a program is used to process and estimate $F$. %, as explained in a previous work~\cite{informatics11020040}. 


\subsection{The sigmoidal logistic function}
The sigmoidal logistic model is a function that can be applied to a time-series~\cite{vikova_transition_2023}, which has different applications such as the representation of epidemic curves. It is represented by the following function:

\begin{eqnarray} \label{sec:FBoltz}
	g \left( x \right) = \frac{H}{1+\mathrm{e}^{A-Bx}},
\end{eqnarray}

where $x$ is the number of days since the first case, $g(x)$ is the expected cumulative number of deaths, $H$ is the height of the wave, $A$ and $B$ are constants related to the pandemic, and $\mathrm{e}$ is a mathematical constant (i.e., an irrational and transcendental number approximately equal to 2.718281828459045).

In this work, a new version of the logistic function is used, which we call the logistic2 function or, simply, the logistic function:

\begin{eqnarray} \label{sec:FBoltz2}
	g \left( x \right) = \frac{H}{1+\mathrm{2}^{A-Bx}},
\end{eqnarray}

where $A$ and $B$ are new pandemic-related constants. Instead of the irrational number $\mathrm{e}$ (which cannot be expressed as a fraction $m/n$, where $m, n \in \mathbb{Z}$), we use the natural number 2 as, in computing, we can efficiently calculate $2^k$ if $k \in \mathbb{N}$. 
Even a school student could understand it.

The calculation of high-order derivatives was performed by hand and was checked using the Octave Software (version 9.1.0). By doing so, we achieved better factorization and presentation of the functions.

\subsection{The program}

The program is developed in R (version 4.4.1 - free software environment for statistical computing and graphics), and the integrated development environment (IDE) RStudio Desktop (2024.09.0+375 "Cranberry Hibiscus" release). It is efficient, simple and easy to understand. The algorithm has the following procedures:

\begin{enumerate}
	\item Load libraries (data structures, nonlinear regression, metrics, graphs, neural networks).
	\item Select country and language.
	\item Load data.
	\item Estimate the model (providing starting values).
	\item Show goodness-of-fit reports ($\sigma$ and RMSE).
	\item Export the model $F(x)$, the first derivative $F'(x)$, and the second derivative $F''(x)$.
	\item Generate the graph (configuring color, size, resolution, output file, among others).
\end{enumerate}

%The source code can be freely downloaded from: (they will be provided during review).


\subsection{Statistical analyses} The Pearson product--moment correlation coefficient or Pearson correlation coefficient ($\sigma$ for short) is a measure of linear dependence and its direction between two quantitative variables. If there is a positive correlation, when one variable changes, the other variable changes in the same direction; if there is no correlation (zero), there is no relationship between the variables; and, if there is a negative correlation, the opposite of the positive correlation occurs. 

The root mean square error (RMSE) is defined as the measure of the differences between the values that are predicted by a model and those that are actually observed~\cite{Sharma2022}. RMSE is not the only measure for this purpose~\cite{Chicco2021}. %, and it has certain drawbacks~\cite{Chicco2021}.%Neither is better than the other and both are used.

%It should be highlighted that all models in this study have two variables---namely, the predicted values and the observed data---therefore, any of the goodness-of-fit measures can be used.

All models in this study have two variables---namely, the predicted values and the observed data---therefore, any of the goodness-of-fit measures can be used. 

Nonlinear regression was used. %Se utiliz� la regresi�n no lineal.
For statistical analysis and graph construction, R software was used.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Results} \label{sec:Results}

\subsection{Case studies}
Four countries were considered: Italy (Southern and Western Europe), Brazil (South America), the Swiss Confederation (west-central Europe), and Peru (South America).

\begin{table*}[h]
	\begin{tabularx}{\textwidth}{>{\hsize=.2\hsize\arraybackslash}X >{\hsize=.2\hsize\arraybackslash}X >{\hsize=.2\hsize\arraybackslash}X >{\arraybackslash}X >{\hsize=.2\hsize\arraybackslash}X}
		\hline
		\textbf{Country} & \textbf{Population} & \textbf{Number of deaths} & \textbf{Start and end date of the pandemic (number of days)} & \textbf{Percentage of deaths} \\
		\hline
		Italy          & 59,500,579   & 188,322 & From 21 February 2020 to 9 March 2023 (1113 days) & 0.32 \% \\
		Brazil         & 213,196,304  & 333,188 & From 17 March 2020 to 9 March 2023 (1088 days) & 0.16 \% \\
		Switzerland    & 8,638,613    & 14,210  & From 06 March 2020 to 9 March 2023 (1099 days) & 0.16 \% \\
		Peru           & 33,304,756   & 219,539 & From 06 March 2020 to 9 March 2023 (1099 days) & 0.66 \% \\
		\hline
	\end{tabularx}
	\caption{Population of the study countries (year 2020), the number of deaths, start and end dates of the pandemic (number of days), and the percentage of deaths in relation to the population}
	\label{table:Population}
\end{table*}

Table \ref{table:Population} presents the populations of the study countries (year 2020), the number of deaths, the start and end dates of the pandemic (number of days), and the percentage of deaths in relation to the population. The population data were downloaded from Population Pyramids of the World (https://www.populationpyramid.net; accessed on 1 September 2024). The last column is the percentage of deaths in relation to the population $= \frac{Deaths}{Population}\cdot 100$.

Day one corresponds to the first death. In Italy, the pandemic began in February 2020, while it began in March in Brazil, Switzerland, and Peru. The pandemic lasted until 9 March 2023. Specifically, it lasted 1113 days in Italy, 1088 days in Brazil, and 1099 days in Switzerland and Peru.

From Table \ref{table:Population}, it can be seen that the country that had a considerably higher percentage of deaths, according to the total number of inhabitants, was Peru ($0.66\%$), followed by Italy ($0.32\%$), and ending in a tie between Brazil and Switzerland ($0.16\%$). 

These countries differ, for example, in terms of their location, population, and number of deaths from COVID-19, among other aspects. However, in relation to the start and end dates of the pandemic, they do not differ greatly.

\subsection{Artificial neural networks}
\begin{table*}[!t]
	\begin{tabularx}{\textwidth}{>{\arraybackslash}X >{\arraybackslash}X >{\arraybackslash}X >{\arraybackslash}X >{\arraybackslash}X}
		\hline
		\textbf{Country} & \textbf{Number of days} & \textbf{Correlation coefficient $R$} & \textbf{Determination coefficient $R^2$} & \textbf{Root mean square error} \\
		\hline
		Italy        & 1113 & 0.9999236 & 0.9998472  & 0.713694600 \\
		Brazil       & 1088 & 0.9999440 & 0.9998881  & 0.008530093 \\
		Switzerland  & 1099 & 0.9999472 & 0.9998945  & 0.003714717 \\
		Peru         & 1099 & 0.9998413 & 0.9996826  & 0.009235257 \\
		\hline
	\end{tabularx}
	\caption{Measures of the ANN}
	\label{CorrelationNN}
\end{table*}

In machine learning, a neural network (NN), also called an artificial neural network (ANN), is a model composed of connected units or nodes called neurons. The main feature is their ability to automatically learn and extract characteristics from input data~\cite{Dai2024}.

Long short-term memory (LSTM) is a type of recurrent neural network (RNN) architecture designed to address the limitations of traditional RNNs, especially when it comes to learning long-term dependencies. The main innovation of LSTM is the use of gating mechanisms to control the flow of information through the network. This allows them to maintain and update their internal state over long periods~\cite{Mienye2024,Torcasio2024,Turchetti2024}.

For better results, it is convenient to select the most appropriate architecture and configure its parameters. 

R provides libraries and frameworks that simplify the process of building and training neural networks, including LSTM networks. The Keras library provides a high-level interface for building and training neural networks. To the best of the authors' knowledge, in most previous works, the Python programming language was used.

To avoid biases in time measurement, we used the computer: HP 11th Generation Intel(R) Core\texttrademark{} i7, and development language: R programming language for statistical computing and graphics.

The data were scaled and split: 90\% of the data was used for training and the rest was retained for testing. The LSTM was used with the following characteristics: three layers (1 input neuron, 200 neurons in the hidden layer, and 1 output neuron), the ``adam'' optimizer, loss = ``mse'', epochs = 200, and batch\_size = 1.

On one hand, the RMSE of the ANN was quite acceptable, being less than $0.72$ (on scaled data, see Table \ref{CorrelationNN}). On the other hand, the Pearson correlation coefficients ($\sigma$) were also quite good: those for Italy, Brazil, Switzerland, and Peru were 0.9999236, 0.9999440, 0.9999472, and 0.9998413, respectively (usually varying with each execution). In ANNs, training took a considerable running time (i.e., more than 8 minutes), and their performance can be improved by expanding their parameters, at the cost of additional time. %In contrast, the runtime of the proposed model was less than a third of a second.

For the ANN, the $\sigma$ did not differ much from each other; this is because ANNs adapt well to data. However, sigmoidal logistic models require the data to follow a pattern. %, which explains the certain differences observed  between the $\sigma$s. %and clarifies why Switzerland had the lowest $\sigma$. 

In general, for some countries the data fit better; however, all countries were fit fairly well. %, as they strictly presented sigmoidal behaviour. 

%Unlike previous works, we revealed that the logistic model is efficient and optimal in terms of its running time. Furthermore, its derivative functions are simple and useful; for example, they allow for description of the dimensions and inflection points of the waves.

Finally, with the parameter settings of an ANN, there is not much that can be done.
%Finalmente, con la configuraci�n de los par�metros de una ANN, no es mucho lo que se puede hacer.

\subsection{Calculating the number of operations of sigmoidal functions}

\begin{table*}[h]
	\begin{tabularx}{\textwidth}{>{\hsize=1.2\hsize\arraybackslash}X >{\hsize=.8\hsize\arraybackslash}X >{\hsize=.8\hsize\arraybackslash}X >{\arraybackslash}X >{\hsize=.1\hsize\arraybackslash}X}
		\hline
		\textbf{Sigmoidal function} & \textbf{Number of power operations} & \textbf{Number of multiplication or division operations} & \textbf{Number of addition or subtraction operations} \\
		\hline
		Gompertz function:\newline $G(x) = a\cdot e^{-e^{b-c\cdot x}}$ & 2 & 2 & 1 \\
		\hline
		Hyperbolic tangent function:\newline $tanh(x) = U - \frac{H}{1+e^{2(x-A)}}$   & 1 & 2 & 3   \\
		\hline
		Boltzmann function:\newline $F(x) = \frac{H}{1+e^{(C-x)/D}}$   & 1 & 2 & 2   \\
		\hline			
		Logistic function:\newline $F(x)= \frac{H}{1+\mathrm{e}^{A-Bx}}$ & 1 & 2 & 2  \\
		\hline
	\end{tabularx}
	\caption{Number of operations in sigmoidal functions}
	\label{table:Noperations}
\end{table*}

The execution time of a sigmoidal function depends on the type and number of operations. In general, the order of the execution time of the operations, from smallest to largest, is as follows: (1) addition or subtraction, (2) multiplication or division, and (3) exponentiation with a real number in the exponent; that is, addition/subtraction is calculated faster than the other two, and multiplication/division is generally faster than exponentiation.

Then, for evaluation, the functions were grouped into types: those based on the Gompertz function (representing the double exponential functions), based on the hyperbolic tangent function (representing the extended functions of the logistic function), and based on the logistic function. The number of operations in each group was then calculated, as detailed in Table~\ref{table:Noperations}.

The Gompertz function was found to be the least suitable, as it has two exponential operations, two multiplication operations, and one subtraction operation. Meanwhile, the hyperbolic tangent function has an additional subtraction operation in relation to the logistic function. The Boltzmann function has two divisions (generally, in computing, the division operation costs more than multiplication). Finally, the logistic function was selected as it has the fewest operations. It could barely be equaled by some extended sigmoidal functions and, therefore, it was considered to be optimal in terms of running time.


\subsection{The model}

First of all, it was necessary to select an optimal and efficient sigmoidal function. 
Subsequently, a general model was developed. %to represent all waves of Covid-19. Finally, using non-linear regression, it was applied to four case studies.
%se desarroll� un modelo general para representar todas las olas del Covid-19. Finalmente, utilizando la regresi�n no lineal se aplic� a cuatro casos de estudio.

The fundamental element of the model is the logistic sigmoidal function, equation (\ref{sec:FBoltz2}). Specifically, the model $F(x)$ is a series of ordered elements, that is, a summation of logistic functions, starting with the function $g(x)_1$ for the first wave, $g(x)_2$ for the second wave, and so on up to $g(x)_n$ for the last wave. The basic structure is:

\begin{eqnarray} \label{sec:FFinal0}
	F \left( x \right) = g(x)_1 + g(x)_2 + \cdots + g(x)_i + \cdots + g(x)_n
\end{eqnarray}


It has the following characteristics or advantages: First, it does not require any form of integration or union of functions, and second, it does not need to add a dummy variable to the data or to the function as in previous publications~\cite{Vilca-Huayta2022,informatics11020040}. 

\begin{theorem} The coupled function with $n$ sigmoidal logistic functions is:
	
	\begin{eqnarray} \label{sec:FFinal1}
		F \left( x \right) &=& \frac{H_1}{1+\mathrm{2}^{A_1 - B_1 x}} + \frac{H_2}{1+\mathrm{2}^{A_2 - B_2 x}} + \cdots +\nonumber\\
		&& \frac{H_i}{1+\mathrm{2}^{A_i - B_i x}} + \cdots + \frac{H_n}{1+\mathrm{2}^{A_n - B_n x}} ,
	\end{eqnarray}
	
	where $F(x)$ is the expected cumulative number of deaths, $x$ is the number of days since the first case, $n$ is the number of waves, $1 \leq i\leq n$, $H_i$, is the height of the wave $i$, $A_i$ and $B_i$ are constants related to the pandemic.
	
	Furthermore, $x_0 < x_1 < x_2 < \cdots < x_n$ are consecutive real numbers that determine the intervals of the logistic functions. The function $g(x)_i$ is defined on the interval $x \in [x_{i-1},x_{i}]$; $\forall x>x_i$, $g(x)_i$ tends to $H_i$, and $\forall x<x_{i-1}$, $g(x)_i$ tends to $0$ (limits of the sigmoidal logistic function).
	%Ademas $x_0 < x_1 < x_2 < \cdots < x_n$ son n�meros reales consecutivos que determinan los int�rvalos de cada funci�n log�stica. La funci�n $g(x)_i$ se define en el int�rvalo $x \in [x_{i-1},x_{i}]$; $\forall x>x_i$, $g(x)_i$ tiende a $H_i$, and $\forall x<x_{i-1}$, $g(x)_i$ tiende a $0$ (l�mites de la funci�n log�stica sigmoidal).
\end{theorem}


\begin{proof}
	Let $x_0 < x_1 < x_2 < \cdots < x_n$ be consecutive real numbers determining the intervals of logistic functions. The function $g(x)_1$ is defined on the initial interval $x \in [x_0,x_1]$, and for $x > x_1$, $g(x)_1$ tends to the constant $H_1$ (right asymptote of a sigmoidal function, $\lim_{x \to \infty } g(x)_1 = H_1$).
	In general, the function $g(x)_i$ is defined on the interval $x \in [x_{i-1},x_{i}]$ (the i-th "wave"); $\forall x>x_i$, $g(x)_i$ tends to $H_i$, and $\forall x<x_{i-1}$, $g(x)_i$ tends to $0$.
	%Sean $x_0 < x_1 < x_2 < \cdots < x_n$ n�meros reales consecutivos que determinan los int�rvalos de cada funci�n log�stica. La funci�n $g(x)_1$ est� definida en el int�rvalo inicial $x \in [x_0,x_1]$, y para $x > x_1$, $g(x)_1$ tiende a la contante $H_1$ (as�ntota derecha de funci�n sigmoidal, es decir, $\lim_{x \to \infty } g(x)_1 = H_1$).
	%En general la funci�n $g(x)_i$ se define en el int�rvalo $x \in [x_{i-1},x_{i}]$ (la i-�sima "ola"); $\forall x>x_i$, $g(x)_i$ tiende a $H_i$, and $\forall x<x_{i-1}$, $g(x)_i$ tiende a $0$.
	
	Proof of $F(x)$ by mathematical induction: Base case (trivial), for two functions; if $x < x_0$, then $F(x) \approx 0$; if $x \in [x_0,x_1]$, then $F(x) \approx g(x)_1$; if $x \in [x_1,x_2]$, then $F(x) \approx H_1 + g(x)_2 \approx g(x)_1 + g(x)_2$, and if $x > x_2$, then $F(x) \approx H_1+H_2$.
	Induction hypothesis: $F(x)\approx \sum _{k=1}^{i-1} \left( \frac{H_k}{1+\mathrm{2}^{A_k - B_k x}} \right ),\ \forall x \in [x_0,x_{i-1}]$.
	Then, $F(x)\approx 0,\ \forall x < x_0$; $F(x)\approx \sum _{k=1}^{i} \left( \frac{H_k}{1+\mathrm{2}^{A_k - B_k x}} \right),\ \forall x \in [x_0,x_{i}]$ (by the induction hypothesis); and $F(x)\approx \sum _{k=1}^{i} H_k, \forall x > x_i$, because, $\lim _{x \to -\infty } g(x)_k = 0 , \forall x > x_i$, and $\forall_k > i$.
	
	
	%Demostraci�n de $F(x)$ por inducci�n matem�tica: Caso base (trivial), para dos funciones; if $x < x_0$, then $F(x) \approx 0$; if $x \in [x_0,x_1]$, then $F(x) \approx g(x)_1$; if $x \in [x_1,x_2]$, then $F(x) \approx H_1 + g(x)_2 \approx g(x)_1 + g(x)_2$, and if $x > x_2$, then $F(x) \approx H_1+H_2$.
	%Hip�tesis de inducci�n: $F(x)\approx \sum_{k=1}^{i-1} \left( \frac{H_k}{1+\mathrm{2}^{A_k - B_k x}} \right),\ \forall x \in [x_0,x_{i-1}]$.
	%Luego, $F(x)\approx 0,\ \forall x < x_0$; $F(x)\approx \sum_{k=1}^{i} \left( \frac{H_k}{1+\mathrm{2}^{A_k - B_k x}} \right),\ \forall x \in [x_0,x_{i}]$ (por la hip�tesis de inducci�n); y $F(x)\approx \sum_{k=1}^{i} H_k, \forall x > x_i$, debido a que, $\lim_{x \to -\infty } g(x)_k = 0, \forall x > x_i$, and $\forall_k > i$.
	
	Finally, $F(x)$ is defined by the equation (\ref{sec:FFinal1}).
\end{proof}

The model is understandable and efficient in running time. %And useful for academic purposes as well.

The construction of $F$ extends to other functions, and it is not necessary that each $g(x)$ be defined exclusively on an interval. The only condition is that the limits of $F(x)$ as $x$ tends to $\pm \infty$ are constant.
%La construci�n de $F$ se exitende a otras funciones, y no es necesario que cada $g(x)$ se defina exclusivamente en un int�rvalo. La �nica condici�n es que los l�mites de $F(x)$ cuando x tiende a $\pm \infty$ sean constantes.

\begin{theorem} The coupled function with $n$ functions is: %La funci�n acoplada con $n$ funciones es:
	
	\begin{eqnarray} \label{sec:FFinal1}
		F \left( x \right) = h(x)_1 + h(x)_2 + \cdots + h(x)_i + \cdots + h(x)_n ,
	\end{eqnarray}
	where, each $h(x)_i$ is a function with left and right horizontal asymptotes, that is,
	
	$\lim_{x\to -\infty} h(x) = L$, and $\lim_{x\to \infty} h(x) = R$, where $L$ and $R$ are constants.
	%donde, cada $h(x)_i$ es una funci�n con as�ntotas horizontales izquierda y derecha constantes, es decir,
	%$\lim_{x\to -\infty} h(x) = L$, y $\lim_{x\to \infty} h(x) = R$, donde $L$ y $R$ son constantes.
\end{theorem}

The proof is similar to Theorem 1. If there are two or more functions with total or partial intersection of their intervals, that is, they are defined in the same interval, they are added.
%La demostraci�n es similar al Teorema 1. Si hay dos o m�s funciones con intersecci�n total o parcial de sus intervalos, es decir, que se definen en un mismo int�rvalo, se suman.

Figure \ref{Fig:Coupled} shows the coupling of four functions: Gompertz, Logistic, algebraic, and wave function. Only the first two are properly sigmoidal, the third is an algebraic function with sigmoidal behavior, and the last is the wave function. It is possible to include different functions that meet the conditions explained, for example, the normal distribution function. On the other hand, the quadratic function $y=x^2$, cannot be coupled.
%La Figura \ref{Fig:Coupled} muestra el acoplamiento de cuatro funciones: Gompertz, Logistic,  algebraica, y funci�n de onda. S�lo las dos primeras son sigmoidales propiamente dichas, la tercera es una funci�n algebraica con comportamiento sigmoidal, y la �ltima es la funci�n de onda. Es posible incluir distintas funciones que cumplan las condiciones explicadas, por ejemplo, la funci�n de distribuci�n normal. Por otro lado, la fuci�n cuadr�tica $y=x^2$, no se puede acoplar.

\begin{figure*}[!t]
	\centering
	\includegraphics[width=410px]{zzPlotZ5.jpg} %\textwidth
	\caption{The coupled function $F(x)$ is made up of four functions: Gompertz, Logistic, algebraic, and wave. Only the first two are properly sigmoidal, the third is an algebraic function with sigmoidal behavior and the last is the wave function.}
	\label{Fig:Coupled}
\end{figure*}

The function $y = F(x)$ is described by equation (\ref{EqCoupled}):

\begin{eqnarray}\label{EqCoupled}
	y &=& \underbrace{2\cdot e^{ -e^{2-x} }}_{Gompertz\ f.} + \underbrace{\frac{2}{1+e^{10-x}}}_{Sigmoidal\ f.} +  \underbrace{\frac{(x-20)}{1+\sqrt{1+(x-20)^2 }}}_{Algebraic\ function} \nonumber\\
	&& - \underbrace{ 4 e^{ - \frac{(x-30)^2}{4} } sin \left( \frac{x-30}{4} \right) }_{Wave\ function} + 1
\end{eqnarray}

The wave function satisfies the requirement of Theorem 2. 
Proof: By definition $-1 \leq sin(x) \leq 1$, it follows that, $- \frac{1}{e^{x^2}} \leq \frac{sin(x)}{e^{x^2}} \leq \frac{1}{e^{x^2}}$. Since the limits of both functions are equal, $\lim_{x \to \infty}- \frac{1}{e^{x^2}}= \lim_{x \to \infty}\frac{1}{e^{x^2}}=0$, it follows that $\lim_{x \to \infty}e^{-x^2} sin \left( x \right) = 0$ (squeeze theorem), and in the same way we can obtain $\lim_{x \to -\infty}e^{-x^2} sin \left( x \right)=0$.\qedsymbol

Finally, the code to graph $F(x)$ in R is:
%Finalmente, el c�digo para graficar $F(x)$ en R es:

\begin{lstlisting}[language=R]
 y = function(x){
  2*e^( -e^(2-x) ) + 2/( 1+e^( 10-x ) ) + 
  (x-20)/sqrt( 1+(x-20)^2 ) -
  4*e^(-(x-30)^2/4)*sin((x-30)/4) + 1
 }
 plot(y, 0, 37, col="green")
\end{lstlisting}

\subsection{The Logistic function of the cumulative number of deaths}

The logistic function $F(x)$ is a set of functions, which are ordered sequentially and do not overlap. Each function corresponds to a wave of COVID-19.

Figure \ref{f:Fig1} represents the cumulative number of deaths from COVID-19. The dataset (observed data) are represented in gray, and the logistic function in blue. Remark, while the sigmoidal shape of the later waves cannot be distinguished due to the scale (thousands on the ordinate axis), it is visible if these (one or two) waves are plotted exclusively.

\begin{figure*}[!t]
\begin{center}
	\subfloat{\includegraphics[width=380pt]{zzPlot9Italy.jpg}
		\label{fig_first_case}} 
	\hfil
	\subfloat{\includegraphics[width=380pt]{zzPlot9Brazil.jpg}
		\label{fig_3_case}}
	\hfil
	\subfloat{\includegraphics[width=380pt]{zzPlot9Switzerland.jpg}
		\label{fig_4_case}}
	\hfil
	\subfloat{\includegraphics[width=380pt]{zzPlot9Peru.jpg}
		\label{fig_5_case}}
	\caption{Cumulative number of deaths from COVID-19 (in thousands) by day.}
	\label{f:Fig1}
\end{center}	
\end{figure*}

\begin{table*}[h]
	\begin{tabularx}{\textwidth}{>{\arraybackslash}X >{\arraybackslash}X >{\arraybackslash}X >{\arraybackslash}X >{\arraybackslash}X}
		\hline
		\textbf{Country} & \textbf{Number of days} & \textbf{Pearson Correlation coefficient $\sigma$} & \textbf{Coefficient of Determination} & \textbf{Root mean square error} \\
		\hline
		Italy        & 1113 & 0.9994238 & 0.9988479  & 1.989923 \\
		Brazil       & 1088 & 0.9999283 & 0.9998566  & 3.013632 \\
		Switzerland  & 1099 & 0.9999178 & 0.9998356  & 0.061266 \\	
		Peru         & 1099 & 0.9999050 & 0.99981    & 1.022324 \\
		\hline
	\end{tabularx}
	\caption{Function adjustment measures.}
	\label{table:Correlation}
\end{table*}

In all cases, the first two or three waves were relatively larger, and they lasted longer (their length on the x-axis is the duration in days); in addition, they had a greater number of victims (represented by the height on the ordinate axis). Meanwhile, subsequent waves were considerably shorter in duration and had fewer victims. In each country, the last wave ends almost flat (the slope tends to zero); under these circumstances, the pandemic ended.

The Pearson correlation coefficients were greater than $R=0.999$ (see the third column of Table~\ref{table:Correlation}), indicating very strong positive linear associations; that is, a high degree of association between the observed data and the function. In the scatter plot, the observed data are very close to the curve. The best adjustment was obtained for Brazil while that for Italy was the worst, but the differences were minimal. The coefficients of determination (fourth column) indicate that the functions obtained are explained or determined by the differences in the days.



The results are very good. Finally, the following hypothesis was accepted: The existence of a correlation that is greater than zero or a statistically significant correlation (with p-value $= 2.2\times10^{-16}$, less than the significance level of $0.01$). %Note: Changing the base $2$ to $e$ or another does not change the correlation or determination coefficients.

Equations (\ref{EqF3}), (\ref{EqF2}), (\ref{EqF4}), and (\ref{EqF5}) represent the functions $F(x)$:


\begin{eqnarray}\label{EqF3}
	\underbrace{F(x)}_{Italy} &=& \frac{31.64}{1+2^{6.6-0.141*x }} + \frac{98.11}{1+2^{10.94-0.033*x }} + \nonumber\\ &&\frac{35.82}{1+2^{32.31-0.045*x }}  + \frac{12.5}{1+2^{45.46-0.051*x }}  + \nonumber\\
	&& \frac{10.71}{1+2^{57.88-0.056*x }} 
\end{eqnarray}


\begin{eqnarray}\label{EqF2}
	\underbrace{F(x)}_{Brazil} &=&  \frac{152.06}{1+2^{5.62-0.045*x }}  +  \frac{464.04}{1+2^{11.9-0.029*x }}  +\nonumber \\ &&\frac{48.72}{1+2^{51.47-0.074*x }}  + \frac{21.45}{1+2^{56.44-0.066*x }}  + \nonumber\\
	&&\frac{14.2}{1+2^{50.35-0.049*x }} 
\end{eqnarray}


\begin{eqnarray}\label{EqF4}
	\underbrace{F(x)}_{Switzerland} &=&  \frac{1.77}{1+2^{6.49-0.192*x }}  +  \frac{7.51}{1+2^{22.68-0.082*x }}  + \nonumber\\
	&& \frac{1.5}{1+2^{14.71-0.041*x }}  + \frac{3.1}{1+2^{24.01-0.036*x }}  + \nonumber\\
	&& \frac{0.36}{1+2^{33.8-0.035*x }} 
\end{eqnarray}


\begin{eqnarray}\label{EqF5}
	\underbrace{F(x)}_{Peru} &=&  \frac{85.94}{1+2^{6.07-0.048*x }}  +  \frac{114.1}{1+2^{15.96-0.041*x }}  + \nonumber\\
	&& \frac{13.22}{1+2^{44.35-0.063*x }}  + \frac{3.4}{1+2^{73.16-0.082*x }}  + \nonumber\\
	&& \frac{3.33}{1+2^{46.8-0.045*x }} 
\end{eqnarray}


%%%%%%%%%%%%%%%%%
\subsection{First derivative}
\label{subsection:1der}

The previous subsection was limited to describing the waves through visual means. To better interpret the mortality rate---represented by the slope or gradient of the curves---the first and second derivative of the functions $F$ were calculated. The third derivative was also calculated, in order to confirm or check the inflection points.
In this way, a correct and elegant interpretation of the results was possible \cite{heusser_baroreflex_2021}.

The first derivative of the Logistic function $F$ or equation (\ref{sec:FBoltz2}) is given in the equation (\ref{sec:FD1}):

\begin{eqnarray} \label{sec:FD1}
	F ' \left( x \right) = \frac{B\,H\,\ln\left(2\right) \, 2^{A-B\,x} }{{\left(1+2^{A-B\,x}\right)}^2},
\end{eqnarray}
where $x$ is the number of days since the first case, $F'(x)$ is the death rate, $H$ is the height of the wave, $A$ and $B$ are the pandemic constants, and $\ln$ is the natural logarithm of base $e$; that is, $\ln(x)=\log_e(x)$.

The first derivative is useful as it represents the slope of the curve or the speed of spread. It is also used to calculate the second derivative. 




%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Second derivative}
\label{subsection:2der}

To obtain the inflection point, it is necessary to calculate the second derivative.

The equation (\ref{sec:FD2}) is the second derivative of the logistic function $F(x)$; that is, the derivative of the equation (\ref{sec:FD1}):

\begin{eqnarray} \label{sec:FD2}
	F'' \left( x \right) = -\frac{B^2\,H\,{\ln\left(2\right)}^2 \, 2^{A-Bx}\,\left(1-2^{A-Bx} \right)}{{\left(1+2^{A-Bx}\right)}^3},
\end{eqnarray}
where $x$ is the number of days since the first case and $F''(x)$ is the second derivative (also called the acceleration) of the cumulative number of deaths.

\subsection{Third derivative}
\label{subsection:3der}

The third derivative $F''' \left( x \right)$ is:

\begin{eqnarray} \label{sec:FD3}
	%F''' \left( x \right) = \\
	\frac{ B^3\,H\,{\ln\left(2\right)}^3\, 2^{A-Bx}\, \left(1-4 \cdot \,2^{A-Bx}+2^{2\,(A-Bx)}\right)}{{\left(1+2^{A-Bx}\right)}^4}
\end{eqnarray}


\subsection{Inflection point}
\label{subsection:pInflection}

A point $c$ on a curve $y = F(x)$ is called an inflection point if $F$ is continuous there and the curve changes from concave upward to concave downward (or vice versa) at $c$~\cite{gu_inflection_2021,hong_enhanced_2022}.

To obtain the inflection point(s), it is first necessary to find the critical point of $F'(x)$; that is, to solve the following equation:

\begin{eqnarray} \label{sec:CPI}
	F''(x) = 0.
\end{eqnarray}

The solution to the equation (\ref{sec:CPI}) is given by $x_0=\frac{A}{B}$ and $y_0=\frac{H}{2}$. 

Finally, $F'''(x_0)$ is evaluated as follows:

\begin{eqnarray} \label{sec:ITest}
	F'''(x_0) = -\frac{B^3\,H\,{\ln\left(2\right)}^3}{8}, \nonumber \\
	F'''(x_0) = - 0.0416 \, B^3\,H.
\end{eqnarray}

The result is different from zero (with $B\neq 0$ and $H\neq 0$), which indicates that $(\frac{A}{B},\frac{H}{2})$ is an inflection point.  It should be highlighted that, according to the method of Higher Derivatives, if $F''(x_0)=0$ holds, then $x_0$ must be substituted into the third derivative. Then, if $F'''(x_0) \neq 0$ holds, then it is an inflection point. %pag. 444

\subsection{Applying high order derivatives}

We continued by deriving the functions $F(x)$. The results are equations (\ref{EqD3}), (\ref{EqD2}), (\ref{EqD4}), and (\ref{EqD5}), represent the first derivatives for Italy, Mexico, Brazil,  Switzerland, and Peru, respectively, using  equation (\ref{sec:FD1}).


	\begin{eqnarray}\label{EqD3}
		\underbrace{F'(x)}_{Italy} &=& \frac{3.1\cdot 2^{6.6-0.141*x}}{(1+2^{6.6-0.141*x})^2 } + \frac{2.22\cdot 2^{10.94-0.033*x}}{(1+2^{10.94-0.033*x})^2 } \nonumber\\
		&& +\frac{1.12\cdot 2^{32.31-0.045*x}}{(1+2^{32.31-0.045*x})^2 } + \\
		&& \frac{0.44\cdot 2^{45.46-0.051*x}}{(1+2^{45.46-0.051*x})^2 } 
		+\frac{0.42\cdot 2^{57.88-0.056*x}}{(1+2^{57.88-0.056*x})^2 } \nonumber
	\end{eqnarray}


	\begin{eqnarray}\label{EqD2}
		\underbrace{F'(x)}_{Brazil}&=& \frac{4.76\cdot 2^{5.62-0.045*x}}{(1+2^{5.62-0.045*x})^2 } + 
		\frac{9.44\cdot 2^{11.9-0.029*x}}{(1+2^{11.9-0.029*x})^2 } \nonumber\\
		&& + \frac{2.48\cdot 2^{51.47-0.074*x}}{(1+2^{51.47-0.074*x})^2 } + \\
		&&\frac{0.98\cdot 2^{56.44-0.066*x}}{(1+2^{56.44-0.066*x})^2 } + \frac{0.49\cdot 2^{50.35-0.049*x}}{(1+2^{50.35-0.049*x})^2 } \nonumber
	\end{eqnarray}


	\begin{eqnarray}\label{EqD4}
		\underbrace{F'(x)}_{ \substack{\text{Swit-} \\ \text{zerland}}}&=& \frac{0.24\cdot 2^{6.49-0.192*x}}{(1+2^{6.49-0.192*x})^2 } +  \frac{0.43\cdot 2^{22.68-0.082*x}}{(1+2^{22.68-0.082*x})^2 } \nonumber\\
		&&+ \frac{0.04\cdot 2^{14.71-0.041*x}}{(1+2^{14.71-0.041*x})^2 } + \\
		&& \frac{0.08\cdot 2^{24.01-0.036*x}}{(1+2^{24.01-0.036*x})^2 } + \frac{0.01\cdot 2^{33.8-0.035*x}}{(1+2^{33.8-0.035*x})^2 } \nonumber
	\end{eqnarray}


	\begin{eqnarray}\label{EqD5}
		\underbrace{F'(x)}_{Peru} &=& \frac{2.87\cdot 2^{6.07-0.048*x}}{(1+2^{6.07-0.048*x})^2 } +  \frac{3.24\cdot 2^{15.96-0.041*x}}{(1+2^{15.96-0.041*x})^2 } \nonumber\\
		&& + \frac{0.58\cdot 2^{44.35-0.063*x}}{(1+2^{44.35-0.063*x})^2 } + \\
		&& \frac{0.19\cdot 2^{73.16-0.082*x}}{(1+2^{73.16-0.082*x})^2 } + \frac{0.1\cdot 2^{46.8-0.045*x}}{(1+2^{46.8-0.045*x})^2 } \nonumber
	\end{eqnarray}


To calculate the second derivative of the function $F(x)$, the equation (\ref{sec:FD2}) was used.


\begin{table*}[h]
	\centering
\begin{tabularx}{\textwidth}{>{\arraybackslash}X >{\arraybackslash}X >{\arraybackslash}X >{\arraybackslash}X >{\arraybackslash}X}
		\hline
		\textbf{Wave}  & \textbf{Brazil} & \textbf{Italy} & \textbf{Switzerland} & \textbf{Peru} \\
		\hline
		1  & 124.38  & 46.68   & 33.87  & 126.08\\%
		2  & 405.60  & 335.82  & 277.2  & 389.42\\%
		3  & 699.89  & 716.04  & 358.53 & 699.39\\%
		4  & 854.16  & 888.31  & 676.27 & 892.89\\%
		5  & 1018.16 & 1034.38 & 964.09 & 1037.25\\%
		\hline
	\end{tabularx}
	\caption{Inflection points by country.}
	\label{table:TableInfle}
\end{table*}


The inflection points play an important role, as they can be considered as so-called ``happy days;'' that is, the day from which the infection rate begins to decrease (mathematically, the slope that was growing begins to decrease). On the graph, this is the point where it changes from concave up to concave down (or vice versa, in other circumstances).

Finally, the inflection points are detailed in Table \ref{table:TableInfle}.

\subsection{Interpretation}

The first derivative of the function $F(x)$ allows us to describe the behavior of the rate or speed of the number of deaths, and is also understood as the slope or gradient of the function $F(x)$.

Figure~\ref{Fig:D1} shows the graphs of the first derivatives of $F(x)$ for each country as blue curves, where the inflection points are indicated by red vertical lines. For comparison between waves or countries, both axes have the same scale. It can be observed that each wave begins with an increasing rate, reaches an inflection point (change), and then begins to decrease. This provides an elegant and analytically based justification for the term "wave".


\begin{definition}[$x$]
	The value of the abscissa $x$ are the days elapsed since the first case or cases. Then, the range (i.e., the number of days from an initial value $x_i$ to another final value $x_j$, with $x_i \leq x_j$) is the duration or the number of days ($x_j - x_i +1$). A full wave, $x$ is the width (or the duration in days). This is valid with respect to the function $F(x)$, the first derivative $F'(x)$, all of the derivatives in general, and their respective graphs.
\end{definition}

\begin{definition}[$F'(x)$]
	The value of the ordinate $F'(x)$ is the slope of the original curve $F(x)$ at point $x$.
	In the case study, the slope is the change of the death rate at the point $x$ (i.e., on day $x$, if $x\in \mathbb{N}$). It is also the velocity and, so, if this value is higher than that on another day, it can be immediately understood as indicating higher lethality. 
\end{definition}

\begin{definition}[$F''(x)$]
	The value of the ordinate $F''(x)$ is the change of rate of $F'(x)$ at point $x$. In the case study, is the acceleration at the point $x$ (i.e., on day $x$, if $x\in \mathbb{N}$). 
\end{definition}

\begin{figure*}[!t]
	\centering
	\includegraphics[width=380px]{zz1DePlotItaly1.jpg} %\textwidth
	%\includegraphics[width=380px]{zzIjaPlotDeMexico1b.jpg}
	\includegraphics[width=380px]{zz1DePlotBrazil1.jpg} %260px
	\includegraphics[width=380px]{zz1DePlotSwitzerland1.jpg} %260px
	\includegraphics[width=380px]{zz1DePlotPeru1.jpg} %260px
	\caption{The first derivative of the COVID-19 logistic function (blue curve) and the vertical projection of the inflection points (red dashed line).}
	\label{Fig:D1}
\end{figure*}

The second wave in Brazil was the largest (in the first figure, it can be observed that it led to more than 400,000 deaths), and its greater height implies that it had a greater lethality.

The first waves in Italy and Switzerland reached a lower death rate and shorter duration, which may be for different reasons, including strict social isolation measures and population differences. In Peru and Brazil, it lasted more than 200 days and had fatality rates below and above one, respectively. Using the proposed model, it was possible to characterize the waves of COVID-19, and it is expected that it can be extended to other countries (which also present similar patterns, according to our pilot study).

Figure~\ref{Fig:D2} represents the graphs of the second derivative of $F(x)$ (for Italy and Brazil only) as blue curves, with the inflection points indicated by dashed vertical lines in red. It can be understood as the acceleration of the number of deaths accumulated due to COVID-19. As the second derivative represents the concavity, when the graph crosses the abscissa axis, it indicates that point as an inflection point; here, the second derivative must change sign, which is interpreted as a change in concavity (positive indicates upward concavity and negative indicates downward concavity). Finally, the second derivatives $F''(x)$ are given as equations (\ref{2DerivativeI}), and (\ref{2DerivativeB}). % (\ref{2DerivativeM}), , and (\ref{2DerivativeP})

\begin{figure*}[!t]
\centering
\includegraphics[width=360px]{zz2DePlotItaly2.jpg}
\begin{comment}
	\includegraphics[width=320px]{zzIjaPlotDeMexico2.jpg}
\end{comment}
\includegraphics[width=360px]{zz2DePlotBrazil2.jpg}
%\includegraphics[width=360px]{zz2DePlotSwitzerland2.jpg}
%\includegraphics[width=360px]{zz2DePlotPeru2.jpg}
\caption{Graph of the second derivative of the logistic function of cumulative COVID-19 deaths in blue (Italy and Brazil) and the projections of the inflection points as red vertical lines.}
\label{Fig:D2}
\end{figure*}


\begin{eqnarray} \label{2DerivativeI}
	\underbrace{F''(x)}_{Italy} &=& 
	\frac{0.3\cdot 2^{6.6-0.141*x}(1-2^{6.6-0.141*x})}{(1+2^{6.6-0.141*x})^3 } + \\
	&&\frac{0.05\cdot 2^{10.94-0.033*x}(1-2^{10.94-0.033*x})}{(1+2^{10.94-0.033*x})^3 } + \nonumber\\
	&&\frac{0.04\cdot 2^{32.31-0.045*x}(1-2^{32.31-0.045*x})}{(1+2^{32.31-0.045*x})^3 } + \nonumber\\
	&&\frac{0.02\cdot 2^{45.46-0.051*x}(1-2^{45.46-0.051*x})}{(1+2^{45.46-0.051*x})^3 } + \nonumber\\
	&&\frac{0.02\cdot 2^{57.88-0.056*x}(1-2^{57.88-0.056*x})}{(1+2^{57.88-0.056*x})^3 } \nonumber
\end{eqnarray}


\begin{eqnarray} \label{2DerivativeB}
	\underbrace{F''(x)}_{Brazil} &= &
	\frac{0.15\cdot 2^{5.62-0.045*x}(1-2^{5.62-0.045*x})}{(1+2^{5.62-0.045*x})^3 } + \\
	&&\frac{0.19\cdot 2^{11.9-0.029*x}(1-2^{11.9-0.029*x})}{(1+2^{11.9-0.029*x})^3 } + \nonumber\\
	&&\frac{0.13\cdot 2^{51.47-0.074*x}(1-2^{51.47-0.074*x})}{(1+2^{51.47-0.074*x})^3 } + \nonumber\\
	&&\frac{0.05\cdot 2^{56.44-0.066*x}(1-2^{56.44-0.066*x})}{(1+2^{56.44-0.066*x})^3 } + \nonumber\\
	&&\frac{0.02\cdot 2^{50.35-0.049*x}(1-2^{50.35-0.049*x})}{(1+2^{50.35-0.049*x})^3 } \nonumber
\end{eqnarray}


%\medskip
%\begin{adjustwidth}{-\extralength}{0cm}
%	\begin{eqnarray} \label{2DerivativeS}
%		\underbrace{F''(x)}_{Switzerland} &= &
%		\left\{ \frac{0.03\cdot 2^{6.49-0.192*x}(1-2^{6.49-0.192*x})}{(1+2^{6.49-0.192*x})^3 } \right\} + \left\{ \frac{0.02\cdot 2^{22.68-0.082*x}(1-2^{22.68-0.082*x})}{(1+2^{22.68-0.082*x})^3 } \right\} + \nonumber\\
%		&&\left\{ \frac{0\cdot 2^{14.71-0.041*x}(1-2^{14.71-0.041*x})}{(1+2^{14.71-0.041*x})^3 } \right\} + \left\{ \frac{0\cdot 2^{24.01-0.036*x}(1-2^{24.01-0.036*x})}{(1+2^{24.01-0.036*x})^3 } \right\} + \nonumber\\
%		&&\left\{ \frac{0\cdot 2^{33.8-0.035*x}(1-2^{33.8-0.035*x})}{(1+2^{33.8-0.035*x})^3 } \right\}
%	\end{eqnarray}
%\end{adjustwidth}
%\medskip

%\medskip
%\begin{adjustwidth}{-\extralength}{0cm}
%	\begin{eqnarray} \label{2DerivativeP}
%		\underbrace{F''(x)}_{Peru} &= &
%		\frac{0.1\cdot 2^{6.07-0.048*x}(1-2^{6.07-0.048*x})}{(1+2^{6.07-0.048*x})^3 } + \frac{0.09\cdot 2^{15.96-0.041*x}(1-2^{15.96-0.041*x})}{(1+2^{15.96-0.041*x})^3 } + \nonumber\\
%		&&\frac{0.03\cdot 2^{44.35-0.063*x}(1-2^{44.35-0.063*x})}{(1+2^{44.35-0.063*x})^3 } +  \frac{0.01\cdot 2^{73.16-0.082*x}(1-2^{73.16-0.082*x})}{(1+2^{73.16-0.082*x})^3 } + \nonumber\\
%		&&\frac{0\cdot 2^{46.8-0.045*x}(1-2^{46.8-0.045*x})}{(1+2^{46.8-0.045*x})^3 } 
%	\end{eqnarray}
%\end{adjustwidth}
%\medskip



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion}
\label{sec:conclusionAndFutureWork}

In almost all countries, and particularly in the study countries (i.e., Italy, Brazil, Switzerland, and Peru), the COVID-19 waves presented a sigmoidal behavior. Therefore, the proposed model obtained high Pearson correlation coefficients.

Through analysis of the type and quantity of operations, it was shown that the logistic function is optimal; that is, it has the shortest execution time compared to similar functions and it is not possible to reduce its operations. Additionally, a new version of the logistic function $g(x)$ was obtained and used, without the addition of operations.

The first three derivatives of the sigmoidal logistic function were calculated, allowing for a precise and elegant analysis. Moreover, the inflection points were obtained (i.e., the point at which the concavity of the function $F$ changes). The derivatives allowed the waves of COVID-19 to be described, analyzed, and compared, as well as providing justification for the well-known term ``waves.''

Highly representative and expressive sigmoidal logistic functions were obtained, which stand out for their efficiency in running time, far surpassing that of the tested artificial neural network. A good fit was obtained with the function, as the waves presented sigmoidal behavior which is typical of epidemics or pandemics.

%On the other hand, the artificial neural network obtained slightly higher Pearson correlation coefficients as, after considerable training, it 
The ANN can adapt to any curve without the requirement that it presents a specific pattern. Therefore, such a model would be preferable if atypical behaviors were present in the observed data, but at the expense of additional time. Also, are suitable for large data sets and high dimensionality in the input variables.

As a final point, neither model is better than the other, and both have their advantages and disadvantages. We consider both to be good for modeling epidemics, particularly when the data follow a sigmoidal pattern.

%Many countries (including the case studies) exhibited five waves, with the first two or three waves being considerably larger than their successors, considering both their duration (in days) and the number of deaths.

The main contribution of this work is the general model that can be used on data that are not necessarily epidemiological; and with any function $h(x)$ with constant left and right asymptotes (that is, $\lim_{x\to -\infty} h(x) = L$, and $\lim_{x\to \infty} h(x) = R$, where $L$ and $R$ are constants, proven by mathematical induction). Also, the method of its construction, and the calculation of high-order derivatives, which allow for the development of an elegant and useful model, as well as justification of the term ``Wave(s)''.

Future work may extend the proposed model to other countries. To build a model $F(x)$, combine different functions (not necessarily sigmoidal). The study of the effects of the waves on different variables, such as the resources used during social isolation measures, is also pending.

\subsection*{Funding statement}
This research was funded by Universidad Nacional del Altiplano (UNA), IPN SIP 20231387 project (RR. N. 1589-2023-R-UNA)

\subsection*{Conflict of interest}
The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

\subsection*{Acknowledgement}
To the administrative staff and the Vice-Rectorate of Research of the Universidad Nacional del Altiplano for their combined efforts and support.


% bibliography style is not predefined; authors can apply their own,
% but it should be sent to matjaz.gams@ijs.si with the accepted paper
%\begin{thebibliography}{99}
%  \bibitem{} Author (year) {\it Title of the book}, Publisher.
%  \bibitem{} Author (year) Title of the paper, {\it Title of the journal}, Publisher, pp. nn--mm.
%  \bibitem{} Author (year) Title of the paper, {\it Title of the proceedings}, Publisher, Location, pp. nn--mm.
%\end{thebibliography}

%\bibliographystyle{plain}
%\bibliographystyle{apalike}
\bibliographystyle{IEEEtran}
\bibliography{IEEEabrv,Paper_format}

\end{document}
