% SIMIODE-TeacherStudent-Template.tex and simiode.cls for LaTeX2e
% Version 22 November 2016.
% 
% This template is to be used with the  class file "simiode.cls". 
%  
% The latest version of this class and template files 
% can be found at
% 
%       <http://www.SIMIODE.org> in the Author section.
%
% A "Quick Start Guide" appears below.
%
% This file must be processed with version LaTeX2e or higher.
% Appreciation to Mark Frantz, IUPUI; John Thoo, Yuba College; and Eric Sullivan, Carroll College.
%
\documentclass{simiode}
%
% With this version of the tex file and the associated SIMIODE.cls file you can produce BOTH a Student Version and a Teacher Version of your modeling scenario or technique narrative.
%
% Student Version contains the title STUDENT VERSION along with the NAME OF SCENARIO, with no abstract, no keywords, no tags, no author identifiers, and no COMMENTS section. The Student Version contains only the STATEMENT of the modeling scenario only with appropriate figures, data, instructions, and perhaps a bibliography appropriate for student resourcing, but nothing in the bibliography which would "give away" the nature of the modeling scenario, etc.
%
%The Teacher Version contains the material in the Student Version, namely the title NAME OF SCENARIO, now headed with   TEACHER VERSION, but also full author  identification, abstract, key words, tag words, and, following the STATEMENT, a section of COMMENTS which can be discussion of presentation  ideas, "solutions," techniques, etc.
%
% Authors can prepare a full Teacher Version using the notions immediately below  to create a pdf for review, but in review please eliminate all information about the author as we use a double-blind referee system. Such information will be restored when accepted and published in SIMIODE at www.simiode.org.
%
%***********   Toggle the following line to turn the teacher version on (1) or off (0). ***********
%
\def\TeacherEdition{0} % Toggle this line's value.
\newcommand{\teacher}[1]{\ifnum\TeacherEdition=1 #1 \fi}
%
% The \teacher{ } command wraps any text that will appear only in the teacher version of
% the document. See this file for examples. This has been done for you in the document below.   
% Rename the pdf after you compile each final version.
%
\usepackage{amsmath,amssymb}
\usepackage{color}
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{url}
%
\renewcommand{\baselinestretch}{1.3}
%
% NOTE regarding macros %%
%
% Any macros defined for your paper should be contained in 
% the top matter. Likewise, any environment definitions such 
% as \newtheorem or \newenvironment should also go in the 
% top matter. 
% 
% Do not \input your macros as separate files within your 
% final version because more files make it harder for the 
% editors and users to keep track of and/or modify and 
% append your paper.
%
% Top matter is above and the beginning of the document below.
%
\begin{document}

\newpage
%
% You only need to change the NAME OF SCENARIO in the next line and delete the footnote. 
% See above for instructions on toggling from Teacher Version to Student Version.
%
\title{\ifnum\TeacherEdition=1 TEACHER VERSION  \else STUDENT VERSION \fi\\World Population }




\markboth{Modeling World Population with Varying Growth Rates}{Modeling World Population with Varying Growth Rates}

\author{Lenka P\v{r}ibylov\'{a},  Jan \v{S}ev\v{c}\'{i}k, Pavel Morcinek and Brian Winkel\\
Masaryk University\\
Brno CZECH REPUBLIC\\
pribylova@math.muni.cz
}

\makeStitlePDFLaTex

\begin{abstract}
We build models of world population using data to estimate growth rate.
\end{abstract}


\section*{SCENARIO DESCRIPTION}

We seek to derive a general model for world population growth using a differential equation.

Let us denote $N(t)$ the size of the world (human)  population at time $t$ in years. Next, let $r(t)$ represent the growth rate of this population (percentage increase of population) at time$ t$. From this we can see the change in the world population  over a period of time $\Delta t$.
\begin{equation}
N(t+\Delta t) = N(t) + r(t)\cdot N(t)\cdot  \Delta t\, .
\label{eq:first}
 \end{equation}

\noindent  {\bf Question 1.} Defend, in words, this model in (\ref{eq:first}).
 
We adjust this relationship to form
\begin{equation}\label{eq:second}
\frac{N(t+\Delta t) - N(t)}{\Delta t} =   r(t)\cdot N(t)\, .
 \end{equation}
and by the limit transition, $\Delta t  \rightarrow 0$, we get the resulting differential equation
 \begin{equation}\label{eq:third}
\frac{dN}{dt} =   r(t)\cdot N(t)\, .
 \end{equation}
 
To simplify the notation, we neglect the time argument of the function $N(t)$ and write $N$. The growth rate $r(t)$ is an unknown function in the variable $t$, which we will estimate using  linear and non-linear regression. 

The data for varying growth rate for world population are shown in Table~\ref{table:RateData}  and displayed in Figure~\ref{fig:RatePlot}. The source is \cite{Meter2020}.


\begin{table}[h]
\begin{center}
\small
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|}\hline
$r_0$  \hspace{-.2 true cm} & \hspace{-.2 true cm}    0.0207	  \hspace{-.2 true cm} & \hspace{-.2 true cm}    0.0197	 \hspace{-.2 true cm} & \hspace{-.2 true cm}  0.0179	 \hspace{-.2 true cm} & \hspace{-.2 true cm}  0.0180	 \hspace{-.2 true cm} & \hspace{-.2 true cm}  0.0181	 \hspace{-.2 true cm} & \hspace{-.2 true cm}  0.0153\\ \hline
World pop  \hspace{-.15 true cm} & \hspace{-.15 true cm}3,700,437,046 \hspace{-.15 true cm} & \hspace{-.15 true cm}4,079,480,606  \hspace{-.15 true cm} & \hspace{-.15 true cm}4,458,003,514  \hspace{-.15 true cm} & \hspace{-.15 true cm}4,870,921,740 \hspace{-.15 true cm} & \hspace{-.15 true cm}5,327,231,061 \hspace{-.15 true cm} & \hspace{-.15 true cm}5,744,212,979\\ \hline
Year  \hspace{-.2 true cm} & \hspace{-.2 true cm}  1970 \hspace{-.2 true cm} & \hspace{-.2 true cm}  	1975 \hspace{-.2 true cm} & \hspace{-.2 true cm}  	1980	 \hspace{-.2 true cm} & \hspace{-.2 true cm}  1985	 \hspace{-.2 true cm} & \hspace{-.2 true cm}  1990	 \hspace{-.2 true cm} & \hspace{-.2 true cm}  1995\\ \hline
	\end{tabular}	
\vspace*{.2 true cm}

\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|}\hline
$r_0$   \hspace{-.2 true cm} & \hspace{-.2 true cm}  0.0133 \hspace{-.2 true cm} & \hspace{-.2 true cm}  	0.0126	 \hspace{-.2 true cm} & \hspace{-.2 true cm}  0.0124 \hspace{-.2 true cm} & \hspace{-.2 true cm}  	0.0119	 \hspace{-.2 true cm} & \hspace{-.2 true cm}  0.0105\\ \hline
World pop \hspace{-.15 true cm} & \hspace{-.15 true cm}6,143,493,823 \hspace{-.15 true cm} & \hspace{-.15 true cm}6,541,907,027 \hspace{-.15 true cm} & \hspace{-.15 true cm}6,956,823,603 \hspace{-.15 true cm} & \hspace{-.15 true cm}7,379,797,139 \hspace{-.15 true cm} & \hspace{-.15 true cm}7,794,798,739\\ \hline
Year 	 \hspace{-.2 true cm} & \hspace{-.2 true cm}  2000	 \hspace{-.2 true cm} & \hspace{-.2 true cm}  2005	 \hspace{-.2 true cm} & \hspace{-.2 true cm}  2010	 \hspace{-.2 true cm} & \hspace{-.2 true cm}  2015 \hspace{-.2 true cm} & \hspace{-.2 true cm}  	2020\\ \hline
	\end{tabular}	
	
	
\caption{\label{table:RateData} Population and varying growth rate for world population over time. Source \cite{Meter2020}.}									
\end{center} 
\end{table}

\begin{figure}[h]
	\begin{center}
		\scalebox{.5}{\includegraphics{RatePlot}}
            	\caption{\label{fig:RatePlot} Plot of the changing annual world population growth rates, $r(t)$ vs.~$t$, time in years. }		
	\end{center}	
\end{figure}		

\noindent  {\bf Question 2.} Assume the growth rate for world population is constant, say $r$. From the data in Table~\ref{table:RateData} offer an estimate for $r$. Create a model for the world population for the future. What does the future hold for the world?

If we believe, based on the data, that the rate at which the world's population grows is NOT constant then we attempt to create a function $r(t)$ for the changing growth rate, certainly based in what we see in Figure~\ref{fig:RatePlot}. We consider the following variations. 



\noindent  {\bf Question 3.} Offer a rationale for each of the two candidates for $r(t)$ the changing growth rate for world population.
\begin{itemize}
\item [(a)]  $r(t) = a_0 + a_1 t$,
\item [(b)]  $r(t) = a_0 + a_1 t + a_2 \sin\left(\frac{2 \pi t}{L} + a_3\right) $\, .   
\end{itemize}

\noindent where $L \approx 30$ indicates the average age of the mother at the birth of one of the offspring.


\noindent  {\bf Question 4.} Show the  solution to (\ref{eq:third}) is (\ref{eq:FirstSol}).

\begin{equation}\label{eq:FirstSol}
N(t) = N(t_0)\cdot e^{\left(\int_{t_0}^t\, r(\tau) d\tau\right)}.
\end{equation}
  

We will assume the present year, $t_0 = 2019$, when there are about 7.714 billion inhabitants on Earth, i.e.~ $N(2019) = 7.714$ billion. 

\subsection*{First model for $r(t)$}

First we consider the model,  $r(t) = a_0 + a_1 t$.  

\noindent  {\bf Question 5.} Show that if the  growth  rate data is fitted to  this simple regression line, the population model takes the form as seen in (\ref{eq:LinearShape}).
\begin{equation}\label{eq:LinearShape}
N(t) = N(t_0) \frac{e^{a_0 t + a_1\frac{t^2}{2}}} {e^{a_0 t_0 + a_1\frac{t_0^2}{2}}}\, .
\end{equation}
 
\noindent  {\bf Question 6.} To estimate the   parameters $a_0$ and $a_1$  use a least squares method to obtain values for the two  parameters, $a_0$ and $a_1$. Use these obtained parameter values for $a_0$ and $a_1$  in (\ref{eq:LinearShape}) and then plot your model for the world's population through  2200. 

The resulting world population model should look something like what appears in Figure~\ref{fig:LinearRatePlot}. Does yours look like this? Explain what this plot is showing. Is this reasonable? Why or why not?  Based on the plot of the world's population in Figure~\ref{fig:LinearRatePlot} determine when the population will peak or reach a maximum.

\begin{figure}[h]
	\begin{center}
		\scalebox{.5}{\includegraphics{WorldPopLinearGrowth}}
            	\parbox{8.5 true cm}{\caption{\label{fig:LinearRatePlot} Plot of the world's population using growth rate (a) $r(t) = a_0 + a_1 t$ with  best values of  $a_0$ and $a_1$.}}
		
	\end{center}
	
\end{figure}	


\subsection*{Second model for $r(t)$}
Next we consider the model $r(t) = a_0 + a_1 t + a_2 \sin\left(\frac{2 \pi t}{L} + a_3\right) $
 
 
\noindent  {\bf Question 7.} Show that if the  growth rate data is fitted to  this more complicated function for $r(t)$, the population model takes the form as seen in (\ref{eq:NonLinearShape}).
\begin{equation}\label{eq:NonLinearShape}
N(t) = N(t_0)  e^{\left(-\frac{-4 \pi  a_0 t+8076
   \pi  a_0-2 \pi  a_1 t^2+8152722 \pi 
   a_1+60 a_2 \cos
   \left(a_3+\frac{\pi  t}{15}\right)+15
   \sqrt{2 \left(5+\sqrt{5}\right)} a_2 \sin
   (a_3)+15 \left(\sqrt{5}-1\right) a_2
   \cos (a3)}{4 \pi }\right)}\, .
\end{equation}
 
\noindent  {\bf Question 8.} To estimate the   parameters $a_0$, $a_1$,  $a_2$, and $a_3$   use a least squares method to obtain values for the four  parameters, $a_0$, $a_1$,  $a_2$, and $a_3$. Use these obtained parameter values for $a_0$, $a_1$,  $a_2$, and $a_3$    in (\ref{eq:NonLinearShape}) and then plot your model for the world's population through  2200. 


The resulting world population model should look something like what appearss in Figure~\ref{fig:NonLinearRatePlot}. Does yours look like this? Explain what this plot is showing. Is this reasonable? Why or why not? Based on the plot of the world's population in Figure~\ref{fig:NonLinearRatePlot} determine when the population will peak or reach a maximum.

\begin{figure}[h]
	\begin{center}
		\scalebox{.5}{\includegraphics{WorldPopNonLinearGrowth}}
            	\parbox{8.5 true cm}{\caption{	\label{fig:NonLinearRatePlot} Plot of the world's population using growth rate (b) $r(t) = a_0 + a_1 t + a_2 \sin\left(\frac{2 \pi t}{L} + a_3\right) $ with your best values of  $a_0$ and $a_1$.}}
	\end{center}
	
\end{figure}	


\noindent  {\bf Question 9.} You should conjecture and support several other functions for $r(t)$ and apply a similar analysis as we do here to predict the future world populations. 

\noindent  {\bf Question 10.} Finally, compare the future for world population as depicted by all four models for $r(t)$: constant growth rate, Model (a), Model (b), and your additional model. What do you think about the future of the world's population based on these analyses?


\section*{More Models for Consideration}
We consider two more models for $r(t)$:
\begin{itemize}
\item [(c)] $r(t) = a_0 + a_1\cdot N(t)$, the logistic model,
\item [(d)] $r(t) = a \left( \frac{N(t)}{A}- 1\right)\left(1 - \frac{N(t)}{K}\right)$, Allee effect model.
\end{itemize}

Model (c), the logistic model, is a model in which the growth coefficient, $r(t) = r(N(t)) = a_0 + a_1\cdot N(t)$, is expressed as a function of the population size, $N(t)$. This gives  the world population growth model of the form seen in (\ref{eq:logistic}):

\begin{equation}\label{eq:logistic}
\frac{dN}{dt} =   r(t)\cdot N(t) =  (a_0 + a_1\cdot N(t))\cdot N(t)\, .
\end{equation}


\noindent  {\bf Question 11.} We will have to estimate the function of the growth coefficient $r(t) =  a_0 + a_1\cdot N(t)$ by linear regression using the data in Table~\ref{table:RateData}. Do that and show that your best fit  coefficients are $a_0= 0.0293993$ and $a_1=-0.00245444$.


\vspace*{.3 true cm}
\begin{figure}[h!]
	\begin{center}
		\scalebox{.5}{\includegraphics{07en.png}}
            	\parbox{8.5 true cm}{\caption{	\label{fig:RatePlotVersusN} Plot of the growth rate estimation (c) $r(t) = a_0 + a_1 N(t) $ with the best estimations of parameters $a_0$ a $a_1$.}}
	\end{center}
	\end{figure}
\vspace*{.3 true cm}

\noindent  {\bf Question 12.}  After substituting Model (c) $r(t) = a_0 + a_1\cdot N(t)$ into the model differential equation (\ref{eq:third}) and solving this differential equation for  the function $N(t)$ verify that we obtain:

\begin{equation}\label{eq:LogisticSolve}        
N(t) = \frac{9.75874\times 10^{23}}{ e^{0.0293993
   t}+3.31938\times 10^{25}}\, .
\end{equation}

\vspace*{.3 true cm}
\begin{figure}[h!]
	\begin{center}
		\scalebox{.5}{\includegraphics{08en.png}}
            	\parbox{8.5 true cm}{\caption{	\label{fig:LogisticModel} Plot of the world's population using growth rate (c) $r(t) = a_0 + a_1 N(t)$ with the best values of $a_0$ and $a_1$.}}
	\end{center}
	\end{figure}
\vspace*{.3 true cm}

\noindent  {\bf Question 13.} Demonstrate the fact that the graph of the world's human population should ``level off'' to 11.8 billion.

The model will never be declining because we chose the declining function $r(N(t))$ (straight line) and thus, the larger the population, the slower it will grow. The smaller the growth coefficient, the smaller the population growth - so growth will be slower, but it will never stop altogether. Such a model can be used, for example, to approximate the capacity of the environment, in our case around 11.8 billion people.

As for Model (d), ``Although the concept of Allee effect had no title at the time, it was first described in the 1930's by its namesake, Warder Clyde Allee. Through experimental studies, Allee was able to demonstrate that goldfish grow more rapidly when there are more individuals within the tank.''\cite{WikiAllee}  This led him to conclude that aggregation can improve the survival rate of individuals, and that cooperation may be crucial in the overall evolution of social structure.  The Allee effect model, has $a$ as the intrinsic rate of increase, $K$ as the carrying capacitiy of the environment, and  $A$ as a critical point or Allee threshold.

We read, ``This graph is an example of a strong Allee effect where the population has a negative growth rate for $0 < N(t) < A$ and are eventually driven to extinction, and a positive growth rate for $A < N(t) < K$ (assuming $0 < A < K$). In cases of a weak Allee effect, populations below the threshold are merely hampered in their rates of growth. The occurrences of Allee effect are more easily observed from small populations as populations above the threshold follow a logistic growth. Since human population has never dropped to such a low level, there is no research concerning about influence of Allee effect to the human population, but Allee effect has been observed in many species, especially animals that hunt for prey or defend against predators as a group.''\cite{Rhea2020}

From a mathematical point of view, we need to estimate $r(N)$ as a quadratic function of $N$ which can be done by using data beyond this threshold\footnote{The threshold is the value, at which a parabola has its maximum.}. However, Table~\ref{table:RateData} does not contain enough information for such an estimation, hence it is necessary to use all 70 measurements of world population size from 1951 to 2020. These data can be found in \cite{Population}. Quadratic growth rate $r(N)$ leads to a world population growth model of the form seen in (\ref{eq:allee}).
\begin{equation}\label{eq:allee}
\frac{dN}{dt} =   a \left( \frac{N(t)}{A}- 1\right)\left(1 - \frac{N(t)}{K}\right) \cdot N(t) \, . 
\end{equation}

\noindent {\bf Question 14.} The function of the growth coefficient $r(t) =a_0 + a_1\cdot N(t) + a_2 \cdot (N(t))^2$ can be estimated by linear regression using the data in \cite{Population} or the data in file population.xlsx. Try it and show that the best fit coefficients are $a_0= 0.016645$, $a_1=1.938243 \cdot 10^{-12}$ and $a_2=-3.693209 \cdot 10^{-22}$. Compare your results with those depicted in Figure~\ref{fig:Allee1}.

\noindent {\bf Question 15.} After substituting Model (d) $r(t) =a_0 + a_1\cdot N(t) + a_2 \cdot (N(t))^2$ into the model differential equation~(\ref{eq:third}) use a numerical method to estimate a solution $N(t)$ of the obtained equation with initial condition $N(2019)=7.714$. Verify that your result agrees with Figure~\ref{fig:Allee2}.

\noindent {\bf Question 16.} Demonstrate the fact that the graph of the world's human population should ``level off'' to 9.83 billion when using Model (d). 

\vspace*{.3 true cm}
\begin{figure}[h!]
	\begin{center}
		\scalebox{.6}{\includegraphics{09en.png}} %upravit obrázek
            	\parbox{8.5 true cm}{\caption{\label{fig:Allee1} Plot of the growth rate estimation (d) $r(t) =a_0 + a_1\cdot N(t) + a_2 \cdot (N(t))^2$ with the best estimations of parameters $a_0$, $a_1$ and $a_2$.}}
	\end{center}	
\end{figure}
\vspace*{.3 true cm}

\vspace*{.3 true cm}
\begin{figure}[h!]
	\begin{center}
		\scalebox{.5}{\includegraphics{10en.png}} %upravit obrázek
            	\parbox{8.5 true cm}{\caption{\label{fig:Allee2} Plot of the world's population using growth rate (d) $r(t) =a_0 + a_1\cdot N(t) + a_2 \cdot (N(t))^2$ with the best values of $a_0$, $a_1$ and $a_2$.}}
	\end{center}	
\end{figure}
\vspace*{.3 true cm}

\begin{thebibliography}{98}

\bibitem{Meter2020} Current World Population. \href{https://www.worldometers.info/world-population}{https://www.worldometers.info/world-population}. Accessed 24 January 2020.

\bibitem{Population} World Population by Year. \href{https://www.worldometers.info/world-population/world-population-by-year/}.  Accessed 24 January 2020.

\bibitem{Rhea2020} Project Rhea. 2020. Logistic Models. \href{https://www.projectrhea.org/rhea/index.php/Logistic\_Models}{https://www.projectrhea.org/rhea/index.php/Logistic\_Models}. Accessed 28 January 2020.

\bibitem{WikiAllee} Wikipedia contributors. (2020, January 25). Allee effect. In {\it Wikipedia, The Free Encyclopedia\/}.  \href{https://en.wikipedia.org/w/index.php?title=Allee\_effect\&oldid=937563935}{https://en.wikipedia.org/w/index.php?title=Allee\_effect\&oldid=937563935}. Accessed 28 January 2020. 

\end{thebibliography}

\teacher{
\section*{NOTES FOR TEACHERS}

\subsection*{Resulting prediction for the case (a) $r(t) = a_0 + a_1 t$} 
After substituting the calculated parameters, we obtain the resulting prediction
 \begin{equation}\label{eq:Result1}
N(t) = 7.714 \frac{e^{0.430088 t - 0.000208\frac{t^2}{2}}} {e^{0.430088\cdot 2019-0.000208\frac{2019^2}{2}}}\, .
\end{equation}

The corresponding development of the human population size is shown in Figure~\ref{fig:LinearRatePlot}.  Based on the model (a) $r(t) = a_0 + a_1 t$, we can predict that the human population will reach its maximum around 2067, at approximately 9.875 billion inhabitants, and then start to decline.

\subsection*{Resulting prediction for the case (b) $r(t) = a_0 + a_1 t + a_2 \sin\left(\frac{2 \pi t}{L} + a_3\right) $} 
The estimation of parameters is $a_0 \doteq 0,460475$,  $a_1 \doteq  -0,223 \cdot 10^{-3}$,  $a_2 \doteq  0,818 \cdot 10^{-3}$ and $a_3 \doteq -0,17292$. After substituting the calculated parameters, we obtain the resulting prediction
 \begin{equation}\label{eq:Result2}
N(t) = 3.5893\cdot 10^{-206}
   e^{-0.000111522
   t^2+0.460475 t-0.000671988
   \sin (0.20944 t)-0.00384732
   \cos (0.20944 t)}
\end{equation}

The corresponding development of the human population size is shown in Figure~\ref{fig:NonLinearRatePlot}.  Based on the model (b) $r(t) = a_0 + a_1 t + a_2 \sin\left(\frac{2 \pi t}{L} + a_3\right) $, we can predict that the human population will reach its maximum around 2061, at approximately 9.717 billion inhabitants, and then start to decline.

\subsection*{Resulting prediction for the case (c) $r(t) = a_0 + a_1\cdot N(t)$}
The estimation of parameters is $a_0= 0,0293993$ and $a_1=-0,00245444$. After substituting the calculated parameters, we obtain the resulting prediction
\begin{equation}\label{eq:Result3}
N(t) =\frac{0,00002284080266}{{ 1,748185649\cdot 10^{20}}\,{{\rm e}^{-
 0,02990992151\,t}}+ 0,000001922576988 }
\end{equation}

The corresponding development of the human population size is shown in Figure~\ref{fig:LogisticModel}. Based on the Model (c), we can predict that the human population will be approaching the approximate value 11.88 billion inhabitants.

\subsection*{Resulting prediction for the case (d) $r(t) = a_0 + a_1 N(t) + a_2 N(t)^2$}
The estimation of parameters is $a_0= 0.016645$, $a_1=1.938243 \cdot 10^{-12}$ and $a_2=-3.693209 \cdot 10^{-22}$. After substituting the calculated parameters, the solution can be approximated by Runge-Kutta method, for example.

The corresponding development of the human population size is shown in Figure~\ref{fig:Allee2}. Based on the Model (d) $r(t) = a_0 + a_1 N(t)+ a_2 N(t)^2$, we can predict that the human population will be approaching the approximate value 9.83 billion inhabitants. The evolution of human population size is subject to weak Allee effect with threshold $A=-4.58$ billion. 

We provide a Mathematics notebook, 1-124-T-Mma-WorldPopulation-TeacherVersion, and a pdf version of the file for those who do not read Mathematica files with analysis of these issues.
}
\end{document}
