Hopf bifurcation in an age-structured prey-predator model with Holling Ⅲ response function

1.   Introduction

In today's world, where significant technological advances have taken places, such as the development of reusable launcher rockets, the design of hypersonic aircraft, and the use of artificial intelligence applications in almost every aspect of our daily lives, so many people have died due to an extremely small entity, called the novel coronavirus (COVID-19). Unfortunately, as of July 13, 2021, the number of cases reported worldwide is 188.021 million, and the total number of deaths is 4.055 million. The number of cases and deaths in Turkey is reported to be 5.486 million and 50,278, respectively, as of the same date.

Just like the other coronaviruses, the COVID-19 is also a single-stranded RNA virus. Although some drugs used in the treatment of diseases transmitted by other RNA viruses are used for the treatment of COVID-19, unfortunately, there is currently no drug developed and proven to be effective for COVID-19. The most common symptoms of COVID-19 are fever, dry cough, joint pain, shortness of breath, bluish lips or face, and decreased sense of taste and smell. Like no other disease in human history, for the COVID-19 pandemic, an incredibly large amount of data is collected, the data collected is constantly updated and publicly shared, and an enormous amount of research is carried out. Probably no one would have thought before that such a micro-entity could bring the whole world together, regardless of language, religion, race, or gender. Researchers from all over the world are investigating many different aspects of the COVID-19 pandemic, such as mathematical modeling of its transmission, vaccination and treatment studies, and its economic, psychological, and sociological consequences.

The mathematical methods used in examining the spread of diseases are generally divided into two classes: statistical methods and mechanistic methods [1]. Pure statistical methods use techniques based on artificial intelligence (AI) or regression. It is clear that the models based on pure statistics do not contain disease-specific information (parameters), so they can only be useful for short-term periods. However, mechanistic methods examine the spread of a specific disease using nonlinear population dynamics. Since these models contain disease-specific information, they allow long-term predictions.

The study [2] conducted by Kermack and McKendrick in 1927 is one of the first studies in which mathematical methods were used to explain epidemiological dynamics. The authors developed a three-equation model, today known as the SIR (susceptible-infected-recovered) model, for simulating the transmission of a contagious disease. Today, most of the deterministic models used for modeling of infectious disease transmission are based on the SIR method, that is, consisting of ordinary differential equations (ODEs). For more on the mathematical models used in biological sciences, one can refer to the materials in [3-13]. The following paragraphs summarize some of the numerous studies conducted on mathematical modeling and numerical simulation of the COVID-19 transmission for many different regions around the world.

Ahmad et al. [14] used a fractional-order mathematical model for simulating the transmission of COVID-19 considering the first 67 days of the breakout in Wuhan. Badr et al. [15] generated a social distancing metric using daily mobility data derived from cell phone data. As a main result of the study, they strongly supported the social distancing and stay-at-home directives in the reduction of case growth rates in US countries. Bertozzi et al. [16] investigated three macroscopic models: exponential growth, self-exciting branching process, and the susceptible-infected-resistant (SIR) model. The authors performed the simulations for California, Indiana, and New York. Cherniha and Davydovych [17] employed a two-equation model for simulating the pandemic for China, Austria, and Poland.

Using the data of March 2020, Coşkun et al. [18] studied the role of climatic factors in the spread of COVID-19 in Turkey. They reported that the population density and wind density are significant factors in the transmission of the virus. Bugalia et al. [19] investigated the transmission dynamics of the virus using a 7-equation model for India. Ivorra et al. [20] developed a mathematical model, called the $\theta$-SEIHRD method, for simulating the transmission of the virus. Besides, they provided an estimation of the need of beds in hospitals in China. Medrek and Pastuszak [21] used artificial intelligence (AI) techniques for simulating the virus spread in Poland, France and Spain. Their model included demographic and geographic factors.

Okuonghae and Omame [22] examined the effects of the measures taken against the spread of COVID-19 for Lagos City of Nigeria using a six-equation model. They reported that if the measures taken were followed by at least 55% of the population, the virus would disappear in the population. Feng et al. used a five-equation model and described how pharmaceutical interventions affect the transmission in [23]. Chen et al. [24] used a bats-hosts-reservoir-people (BHRP) network for modeling the virus transmission in Wuhan City. Samui et al. [25] employed a four-equation fractional-order model for India. Tuan et al. [26] used a six-equation fractional-order model for modelling the virus transmission for the entire world population. Çakan [27] developed a four-equation model taking into account the health care capacity and pointed out that the number of hospital beds, intensive care units, and respirators should be increased in the fight against such pandemic diseases.

Yadav and Verma [28] developed a fractional-order model based on the Caputo–Fabrizio derivative for simulating the virus transmission in Wuhan City. Atangana [29] constructed a mathematical model for simulating the spread in Italy taking into account the effects of lock-down. The author proposed new fractal-fractional differential and integral operators and applied them to the introduced mathematical model. Atangana and Atangana [30] emphasized how simple but effective the face masks are against the coronavirus, and simulated the possible side effects numerically. Khan et al. [31] studied the dynamics of COVID-19 considering the stay-at-home directives for China. They developed the mathematical model using ordinary and fractional-order differential equations in the Atangana–Baleanu sense.

Viguerie et al. [32,33] proposed a mathematical model based on partial differential equations (PDEs), called a spatiotemporal model, by interpreting the transmission of the virus in both space and time, similar to those in reaction-diffusion systems. The authors performed the simulations on the geometry of the Italian region Lombardy. Zhao and Feng employed a seven-equation model for the transmission of COVID-19 reporting that properly designed staggered-release policies can do better than simultaneous-release policies in terms of protecting the most vulnerable members of a population, reducing health risks overall, and increasing economic activity in [34]. In reference [35], the authors considered the dynamics of the COVID-19 epidemic in Ukraine, and modeled it using a modified age-structured compartmental SEIR (susceptible, exposed, infected and recovered) framework. They also provided a short-term forecast of further dynamics of COVID-19. Sarkar et al. [36] employed a six-equation model for simulating the transmission of the pandemic in seventeen different provinces of India. Shaikh et al. [37] proposed a fractional-order model for simulating the transmission of the virus in India using the data of 14 March, 2020 to 26 March, 2020.

Naeem et al. studied the Pythagorean $m$-polar fuzzy topology (P$m$FT) approach with a multi-criteria decision analysis tool, called the technique for order of preference by similarity to ideal solution (TOPSIS), in [38]. In reference [39], Saha et al. formulated a compartmental SEIRS model of the pandemic where a separate equation is incorporated to reflect the information which induces behavioral changes. Moradian et al. discussed the importance of bringing the world's scientists together and interdisciplinary collaborations to find effective solutions for controlling the pandemic in [40]. Mohamed and Rezaei [41] underlined the importance of the combination of the computational drug repurposing methods and clinical trials delivers on the drug discovery phenomena. Cheema et al. [42] examined the effects of the social distancing and quarantine measures taken in Wuhan from 13 February, 2020 to 24 March, 2020. They reported that the quarantine measures prevented the transmission of the virus significantly.

Kuhl summarizes what we have learned about the emergence and spread of COVID-19 in ten items in her elegant study [43]: the virus spreads exponentially when controlled, is contagious like other coronaviruses, as long as an effective vaccine is found and applied to the population correctly, it will not leave the humanity, lower reproduction numbers can be obtained, movement restrictions are highly effective in preventing the spread of the virus, there is a delay of about two weeks between imposing the restrictions and observing the effects, many COVID-19 cases are asymptomatic and not reported, there are some differences between the mechanistic/statistical models and actual data, selective reopen may be more effective than voluntary quarantines, and tests are crucial for a safe reopen.

The topic is so dynamic that new developments in vaccine studies are constantly occurring, the effectiveness of these vaccines in combating new emerging virus variants tends to vary, and the closing and opening measures are taken/removed/updated by countries are constantly changing. As a result of all of these developments, it is unsurprising that the mechanistic and statistical models developed to simulate the spread of the virus are constantly changing and updating. For the most up-to-date information on COVID-19 from a variety of perspectives, see reference [44-52] and the materials therein.

In this study, the mathematical model used to simulate the human-to-human spread of COVID-19 for Wuhan city in [53] is analyzed for the case of Turkey, and the related parameters are calibrated accordingly. Rather than analyzing the mathematical model from theoretical perspectives, we focus on the problem computationally and compare the results visually. The mathematical model, which includes the factor of super-spreaders, consists of eight ODEs and has a nonlinear structure. The Galerkin finite element method (GFEM) is used for solving the governing equations, and the Newton–Raphson method is employed for solving the system of nonlinear algebraic equations resulting from the temporal discretization.

In Section 2, the mathematical model for simulating the transmission of COVID-19 is introduced. Besides, some mathematical properties of the model are provided. In Section 3, the mathematical model is tested and compared with the real data. Section 4 studies a sensitivity analysis of the model parameters to the basic reproduction number, $R_{0}$. In Section 5, the findings are discussed in detail, and a comprehensive outlook is given. Turkey's data covering the first 81 days of the pandemic is presented in Appendix A, and the solver developed for performing the simulations is provided in Appendix B.

2.   The mathematical model for COVID-19 transmission

In this section, we calibrate the parameters used in the mathematical model introduced by Ndaïrou et al. in [53,54] for the transmission of COVID-19 with a case study of Wuhan. In this mathematical model, the constant total population, $N$, is classified into eight epidemiological groups: $S(t)$, $E(t)$, $I(t)$, $P(t)$, $A(t)$, $H(t)$, $R(t)$, and $F(t)$. These epidemiological groups correspond to the susceptible $(S)$, exposed $(E)$, symptomatic and infectious $(I)$, super-spreaders $(P)$, asymptomatic but infectious $(A)$, hospitalized $(H)$, recovery $(R)$, and fatality $(F)$ classes, respectively. The transmission network is shown in Figure 1.

The mathematical model is given as follows:

where the parameters are described and explicitly given in Table 1. In this model, the number of deaths at an instant time (day) is given as

Besides, one should notice that the total population, $N$, is constant and defined as

where the number of fatal patients is very small compared to the total population, $F\ll N$. Since their behavior is not predictable, the transmissibility from asymptomatic individuals is ignored in the governing equations. This is, in fact, currently a controversial issue among epidemiologists [53].

It should be also noticed that the exact solutions to the differential equations (2.5), (2.7) and (2.8) can be determined easily through basic integration. Furthermore, the function $S(t)$ can be expressed in terms of the functions of other compartments using Eq (2.10). Therefore, any qualitative analysis of the Eqs (2.1)–(2.8) can be studied by considering the remaining equations of the system. Let us recast the governing equations in the form of

where the vectors $\bf{F}$ and $\bf{V}$ represent the rate of appearance of new infections and the net rate of transfers in the compartments, respectively. These vectors are defined as

The Jacobians associated with matrices $\bf{F}$ and $\bf{V}$ are given as follows:

Notice that the Jacobian $\bf{J_{F}}$ can be given as

since the total population, $N$, and the number of the individuals in the susceptible class, $S$, are very close initially, at $t = 0$. The inverse of the Jacobian $\bf{J_{V}}$ is given as

where the entry $A$ appearing in the matrix is defined as

Thus, the matrix $\bf{J_{V}}^{-1}\bf{J_{F}}$ is given as follows:

where

Then, the basic reproduction number, $R_{0}$, is defined as given in [55,56]:

where $\omega_{i} = \gamma_{a}+\gamma_{i}+\delta_{i}$, $\omega_{p} = \gamma_{a}+\gamma_{i}+\delta_{p}$, $\omega_{h} = \gamma_{r}+\delta_{h}$, and $\rho$ denotes the spectral radius.

In reference [53], the authors performed a local stability analysis of system (2.1)–(2.8) in terms of the reproduction number, $R_{0}$. Furthermore, they conducted a sensitivity analysis of the parameters. A similar analysis is performed for Turkey's data in Section 4.

One should notice that the basic reproduction number, $R_{0}$, is greater than one means that the average number of individuals in the population infected by an individual with the disease is more than one. This indicates that the disease will continue to spread depending on the magnitude of $R_{0}$. Similarly, but on the contrary, the fact that $R_{0}$ is smaller than one means that the disease will disappear in the population over time. For this reason, the correct determination and interpretation of $R_{0}$ are crucial in combating the disease.

3.   Numerical simulations

In this study, we examine the first wave of the COVID-19 outbreak in Turkey because, during the period between July and December 2020, the Turkish authorities reported the number of hospitalized patients rather than the total cases. Therefore, the data we consider covers a period of 82 days, 13 March, 2020 to 2 June, 2020, and is given in Appendix A.

Turkey had a population of approximately $83,150,000$ in the early 2020. In this case, the total population is taken as $N = 83,150,000 / 1,075$ in a similar fashion as the authors followed in [53]. The denominator in this formula is used to account for the impact of travel restrictions imposed during the epidemic in the mathematical model given by Eqs (2.1)–(2.8).

The system of nonlinear ODEs is solved by using the GFEM. Particularly, we use Lagrange elements of degree $p = 3$. The interested reader can refer to [57-59] for more on the finite element methods. All computations are carried out in the FEniCS environment. The details concerning the FEniCS Project are available in [60-62]. The interested reader can also find a sample solver in Appendix B to reproduce the results and/or perform their own simulations.

The initial conditions are as follows: $S(0) = N-24$, $E(0) = 0$, $I(0) = 4$, $P(0) = 20$, $A(0) = 0$, $H(0) = 0$, $R(0) = 0$, and $F(0) = 0$. Here, the epidemiological class $I$ denotes the symptomatic and infectious individuals and it is reported that four cases were detected on March 13, 2020 [63]. Furthermore, based on the study [53], the number of super-spreaders is set as five times the number of symptomatic and infectious individuals. In other words, it is examined how four infected individuals and twenty super-spreaders in the population affect the others over time. All the parameters used in the simulations are given and described in Table 1. These constant parameters have been obtained by calibrating those given in [53,54] for Wuhan to Turkey's data.

In Figures 2 and 3, the results obtained by solving Eqs (2.1)–(2.8) and the actual data are compared for the number of confirmed cases and number of confirmed deaths, respectively. The reproduction number is computed by using Eq (2.18) and found as $R_{0} = 3.7364$. It means that, for each individual with the disease, the average number of people infected by the individual is $3.7364$. In Figures 4 and 5, simulations for various epidemiological groups are given.

4.   Sensitivity analysis

The sensitivity analysis of the parameters given in the previous section is essential in determination of the dynamics of the pandemic. This analysis can be used to measure the relative change caused by a change in any parameter. The normalized forward sensitivity index of a parameter $\Theta$ with respect to a differentiable quantity $R_{0}$ is given as follows [64,65]:

where $\Theta$ is the parameter under consideration. As one can deduce from Eq (4.1), the sensitivity index of a parameter might be a constant or depend on several parameters involved in the system.

It can be seen in Table 2 that the most sensitive parameters to the reproducing number, $R_{0}$, are $\beta$, $\rho_{1}$, $l$, and $\gamma_{r}$. For example, $\Psi_{\beta}^{R_{0}} = 0.99595$ means an increase (or decrease) in $\beta$ by $10\%$ results in an increase (or decrease) in $R_{0}$ by $9.9595\%$. These results are very similar to those obtained for Wuhan city in [53,54].

5.   Conclusions

In this study, we have applied the mathematical model developed in [53] for simulating the spread of COVID-19 in Wuhan city to Turkey's data. This mathematical model examines the total population, $N$, in eight different groups using a system of ODEs. The GFEM is used to solve the system of ODEs. It was seen that there are some differences between the real data and numerical simulations:

● It has been observed that the spread of the virus was successfully controlled at the beginning of the outbreak:

    – The measures to prevent the pandemic seem taken timely and applied quickly by the authorities,

    – The population behaved in accordance with the measures taken.

● The number of deaths was above the curve in the early stages of the pandemic. This result may be caused by the followings:

    – Early times in the pandemic, the lack of knowledge about the treatment may have resulted in a higher number of deaths,

    – Individuals who contracted the virus and developed symptoms may have mistaken their illness for a simple cold and not sought a treatment.

● Right after the peak point, both the number of cases and deaths were above the curve. Some of the many reasons that might cause this:

    – The quarantine measures may not have been followed by the population or not strictly managed in this period,

    – The drop in the number of cases and deaths towards the summer months may have encouraged people to be more reckless.

Although there is a significant difference between real data and numerical results, this is the case for almost all studies in the literature [43]. Also, as mentioned before, this study has solely focused on the first wave of the COVID-19 outbreak in Turkey. The main reason for this decision is that the Turkish authorities considered it more appropriate to report the number of cases in a different way during a certain period after the first wave. This situation created a limitation for the functioning of the mathematical model we use. Future studies are invited to use the data of the subsequent waves to make comparisons.

Some natural extensions of this study can be listed as follows:

● The human-to-human transmission model used can be extended to different transmission networks, e.g., a bat-host-reservoir-people (BHRP) network can be employed.

● The model can be used to simulate the transmission processes of COVID-19 in different populations, for example, for another country or a city.

● The transmissibility from asymptomatic individuals can be included in the mathematical model.

● The COVID-19 had been known to affect individuals over middle age more than others. Speculative and unproven information on the effects of the virus on individuals of different genders, nationalities, and blood types is constantly circulating. It will pave the way for intriguing studies and information by adding demographic variables to the relevant models.

● To obtain more accurate numerical results, the measures that have been taken to prevent the epidemic can be included in the relevant model. Besides, it is possible to examine the effects and protection of vaccines developed for the virus. Also, the effects of environmental factors such as air humidity, air pollution, and wind intensity on the spread of the virus can be investigated.

● Spatio-temporal models, in which both spatial and temporal simulations of the spread of the virus are performed, can be considered.

Acknowledgments

We would like to express our gratitude to the anonymous reviewers for their insightful comments, which helped us improve the manuscript.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendixs

A.   Turkey's Data

The data covering the period of 82 days, March 13–June 2, 2020 [63]:

and

where $L_{C}$ and $L_{D}$ denote the data for the confirmed cases and deaths, respectively.

B.   FEniCS Code

The solver developed in FEniCS 2019.1.0 is given as follows: