key: cord-0733632-dkxcjfrp
authors: Ghosh, Kalpita; Ghosh, Asim Kumar
title: Study of COVID-19 epidemiological evolution in India with a multi-wave SIR model
date: 2022-04-27
journal: Nonlinear Dyn
DOI: 10.1007/s11071-022-07471-x
sha: 6f63f7644fde5f88cb1dc900467aa6eca5775415
doc_id: 733632
cord_uid: dkxcjfrp

The global pandemic due to the outbreak of COVID-19 ravages the whole world for more than two years in which all the countries are suffering a lot since December 2019. In this article characteristics of a multi-wave SIR model have been studied which successfully explains the features of this pandemic waves in India. Origin of the multi-wave pattern in the solution of this model is explained. Stability of this model has been studied by identifying the equilibrium points as well as by finding the eigenvalues of the corresponding Jacobian matrices. In this model, a finite probability of the recovered people for becoming susceptible again is introduced which is found crucial for obtaining the oscillatory solution in other words. Which on the other hand incorporates the effect of new variants, like delta, omicron, etc in addition to the SARS-CoV-2 virus. The set of differential equations has been solved numerically in order to obtain the variation of susceptible, infected and removed populations with time. In this phenomenological study, some specific sets of parameters are chosen in order to explain the nonperiodic variation of infected population which is found necessary to capture the feature of epidemiological wave prevailing in India. The numerical estimations are compared with the actual cases along with the analytic results.

The outbreak of COVID-19 in December 2019, due to the spreading of infectious coronavirus named SARS-CoV-2 finally triggered to a global pandemic in which all the countries are suffering a lot till to date. This highly contagious disease has been transmitted to millions of people globally where a fraction of infected people is succumbed to it eventually. Most of the countries witnessed multiple peaks of epidemiological infections in its evolution which is counted as number of waves. Among the most affected countries, for example, countries in European union (CEU), USA, Russia and Canada experienced five successive epidemiological waves while the third wave is going on in India, Indonesia and Brazil. Which means that five peaks of the epidemiological infections are found in CEU, USA, Russia and Canada, while three distinct peaks are noted in India, Indonesia and Brazil. However, peaks of epidemiological infections in India are sharper than those found in Brazil. So, it is obvious that characteristics of these epidemiological waves found in different countries are not similar. Widths and heights of the peaks are different as well as the separation between them are not equal. It is observed that the same individual gets infected multiple times during this pandemic in all countries. This is due to fact that new variants of SARS-CoV-2 virus emerge in due course.

In order to study the dynamics of the pandemic in a deterministic approach, several models have been proposed based on the classic SIR model introduced by Kermack and Mckendrick in 1927 [1] . This model yields only a single peak for infection in its evolution, so, henceforth it is referred to as the single-wave SIR (SWSIR) model. Most of the previously studied models are mainly formulated to explain the characteristics of the solitary wave of the pandemic. Those are known as SEIR, SIQR and their hybrids which are all derived from the SWSIR model [2] [3] [4] [5] [6] . In those attempts, effects of quarantine, isolation, latent time of infection and other factors are taken into account where several predictions and forecasts are available [7, 8] . Nonetheless, periodic outbreak of disease in terms of SIR, SIQR, and SEIQR models are studied analytically in great details [9] [10] [11] . Evolution of the COVID-19 pandemic for different countries have been studied with the help of those models where origin of differences in its features are explained. Case studies on few countries are available [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] . However, generation of multi-wave pattern in the COVID-19 pandemic evolution in a deterministic approach through a single model is not reported before.

Hence, a multi-wave SIR (MWSIR) model is formulated by modifying the SWSIR model in order to explain the dynamics of epidemiological infections including the origin of multiple waves found in the infection pattern. This is a fact that the same individual has been infected repeatedly during this pandemic period. So in this model, probability of the same individual for getting infected multiple times after curing is taken into account [22] . Time delay between successive infections of the same individual is introduced in the MWSIR model. In this work, periodic outbreak of COVID-19 in India will be studied, and so, the daily new cases in terms of seven-day moving average in India during January 30, 2020 to March 14, 2022 (775 days) are shown in Fig. 1 . It reveals that the duration of those waves are different, so, it is not at all periodic. The approximate duration of first and second waves are T 1 = 370 and T 2 = 300 days, respectively, while the third wave crosses 105 days (T 3 = 105). Hence, separation between adjacent peaks are different, and in addition, width and height of the peaks are unequal.

In Sect. 2 a MWSIR model is considered which gives rise to almost periodic epidemiological oscillation. The stability of this model is studied in Sect. 3 analytically while its periodic feature is studied numerically in Sect. 4 . In order to produce nonperiodic oscillation, this model has been modified further. The epidemiological oscillation produced in this study is compared with the daily new COVID-19 cases in India. A discussion based on those results is available in Sect. 5.

In order to explain the epidemiological waves of infection due to the spreading of COVID-19 in India, a MWSIR model is formulated where the total population is divided into three dynamic sub-populations, those are known as susceptible, infected and removed. Susceptible population is constituted by the individuals those are otherwise healthy but have a probability for having infected at any time in future. The intermediate stage is comprised of infected population where the individual is instantaneously contagious. In the terminal stage individuals are either recovered or succumbed to the disease, however, they jointly referred to as removed population. As the ongoing COVID-19 pandemic exhibits periodic outbreak of the disease globally, a modified version of SWSIR model is necessary to study the dynamics of this pandemic.

Therefore, in order to understand the nature of periodic outbreak, a MWSIR model has been formulated by modifying the primitive version of SIR model as described below. In this model the recovered population has a nonzero probability for getting infected again. As a result, the susceptible population gets enriched with time which is measured in terms of infected population with a time delay, τ . Thus, the MWSIR model is defined by a set of coupled first-order nonlinear differential Eq. (1):

where S(t), I (t), and R(t) are the number of susceptible, infected and removed people at time t. The positive 

Flowchart for the MWSIR model constants μ, α, β and γ are the rates of birth, transmission (or infection) per individual, removed and resusceptibility, respectively. The model incorporates the fact that the recovered people become susceptible again after the mean period τ [22] . However, I (t − τ ) = 0, when t < τ. The SWSIR model will be restored when μ = 0, and γ = 0. The system in this formulation is closed in a sense that the total population, N = S(t) + I (t) + R(t), does not change with time, which means the rate at which individual suffers natural mortality is also given by the parameter μ. In other words, mean lifespan of individual turns out to be 1/μ. In the same way, 1/β can be regarded as the mean recovery time for the infected individual. The flowchart of this model is shown in Fig. 2 . Now it is pertinent to explain how does this MWSIR model yield multiple peaks in its evolution. However, looking back on to the SWSIR model (μ = 0, γ = 0), Eq. (1a) leads to the fact thatṠ < 0 for all t. It means the number of susceptible individuals is always decreasing irrespective of the initial value, S(0). S(t) is monotone and positive which corresponds to lim t→∞ S(t) = S ∞ . On the other hand, Eq. (1c) corresponds toṘ > 0 for all t, which means the number of removed individuals is always increasing. Although R(t) possesses monotonic behaviour irrespective of the initial conditions but it is bounded by total population N , which corresponds to

Equation (1b) reveals that the number of infected individuals increases monotonically as long as β < αS(t) up to some maximum value and decreases to zero since the opposite condition, β > αS(t), is being satisfied beyond a certain time and thereafter. So, the epidemic dies out finally as lim t→∞ I (t) = I ∞ → 0. Thus it yields a solitary broad peak in the evolution for the number of infected individuals. Equation (1b) further predicts when this maximum will be attained as it can be determined by the conditionİ = 0, which corresponds to the relation, S(t) = β/α. Further it can be shown that the value of the maximum can be expressed as [5] ,

On the other hand, for the MWSIR model (μ = 0, γ = 0), different picture is observed since Eq. (1a) could be decomposed along the time domain aṡ

where the two opposite conditions are successively satisfied with the increase of time. To understand this point, let us start from t = 0. As in the beginning

is satisfied, so,Ṡ < 0. Hence, the number of susceptible individuals decreases monotonically up to some minimum value. ObviouslyṠ > 0 after some time since I (t − τ ) > 0 when t > τ, so, Eq. (3b) is fulfilled and as a consequence, S(t) increases afterwards. However, Eq. (3a) reveals that there is a possibility for the condition,Ṡ < 0, to be satisfied again at a certain time when S(t) becomes very large. Thus the conditions, S < 0 andṠ > 0, are found to be valid alternately with the time leading to oscillations in its evolution. As a result, number of infected and removed individuals also show oscillatory behavior maintaining some phase difference between themselves. So eventually the populations altogether exhibit multiple peaks in their respective evolution which is shown Fig. 3 . Positions of the extrema for the infected population could be determined easily by using Eq. (1b), as this one remains the same both for the SWSIR and MWSIR models. Therefore, positions of the extrema for the infected population could be determined using the same relation, S(t) = β/α, also in this case. So, on the other hand, value of α can be found by locating the position of the peak for the infected population by plotting both the variations of S(t) and I (t) in the same scale and this geometrical method will be described later while confirming the value of α. In addition, value of the first peak has been determined using Eq. (2).

The basic reproduction number, R 0 plays an important role in epidemiology, which is defined as average number of new infections per infected individual. So, in this three-compartment system, it is equal to the transmission rate multiplied by the infectious period, R 0 = α/(β + μ) [4, 11] . The value of R 0 is crucial to have an idea about how does the disease flow in the whole population and at the same time it provides clue to control its spreading.

Theorem: Every solution of the MWSIR model along with its initial conditions is a subset in the interval [0,∞) and {S(t), I (t), R(t)} ≥ 0 for all values 0 ≤ t < ∞. In order to prove this, Eq. (1a),

has been expressed as d dt Se

Hence, upon integrating

As S(0) > 0,

Thus,

As a result, S(t) > 0 at any time. Similarly, it can be shown that I (t) > 0, and R(t) > 0, which complete the proof. In order to determine the stable equilibrium points, the rate of change of S(t), I (t) and R(t) with time are made zero. The trivial solution leads to I (t) = 0 at any time, which corresponds to disease-free equilibrium (DFE). Mathematically, this point is expressed as (S * , I * , R * ) tr = (N , 0, 0), where the entire population becomes susceptible but no infected individual. This is true for both t < τ and t ≥ τ . Anyway, this DFE is insignificant as it does not correspond to dynamics of epidemic by any means.

The nontrivial solution leads to I (t) = 0, at finite time, but, I (t → ∞) → 0, which corresponds to endemic equilibrium and that is observed in SWSIR model. However, in this MWSIR model, nontrivial solutions mean I (t) > 0, at any time. So these dynamic equilibrium points appear periodically with time without showing any signature of ending if there is no damping. These equilibrium points are marked by (S * , I * , R * ) ntr = β+μ α , μN β+μ − μ α , β N β+μ − β+μ α + μ α , when t < τ, however, for t ≥ τ , (S * , I * (t), R * (t)) ntr = β+μ α , μN β+μ + γ I (t−τ ) β+μ − μ α , β N β+μ − β+μ α − γ I (t−τ ) β+μ + μ α . Hence, equilibrium points are static when t < τ, but dynamic when t ≥ τ . In order to understand the dynamics close to the equilibria, Jacobian matrix (J ) is constructed which can be expressed at the respective equilibrium points. Jacobian matrix has the form

The eigenvalues of Jacobian matrix (λ 0 , λ ± ) can be given as

For the DFE, (S * , I * , R * ) tr = (N , 0, 0) , the eigenvalues are ⎧ ⎨

As a result, it will behave like a stable point when α < (β + μ)/N , which is similar to the SWSIR model. It is obvious that the eigenvalues will be complex for the general case when αβ > α 2 (I * − S * ) 2 + β 2 2(S * + I * ) .

be oscillatory around these points [4] . The solutions for S(t), I (t), and R(t) in this linearized model could be found in terms of the eigenvalues and eigenvectors of J , however, they are not necessary at this moment.

As the values of S(t), I (t), and R(t) cannot be determined analytically by solving the set of Eqs. (1), they have been solved numerically. The oscillations of S(t), I (t), and R(t) are noted for long time and it is observed that all are oscillating with almost constant amplitude. Four successive waves are shown in Fig. 3 , where it is observed that the mean time period of oscillations, T , is always greater than τ . The reason behind this oscillation can be understood in terms of Eqs. (3) which was explained before. As no damping factor is introduced here, amplitudes of the respective waves for S(t), I (t), and R(t) do not decay with time.

The time period of oscillations for S(t), I (t), and R(t) are obviously the same, although it decreases with the number of cycles. However the difference, T − τ vanishes with the increase of τ . This feature is shown in Fig. 4 where variations of time period T along with The variation of T with α and β for fixed value of τ is shown in Fig. 5 . It clearly indicates that T > τ, for any values of α and β. Moreover, the difference, T −τ steadily decreases with α but increases gradually with β.

Features of the ongoing COVID-19 epidemiological waves observed in India are shown in Fig. 1 . It reveals that the duration of successive waves are decreasing as T 1 > T 2 > T 3 . On the other hand, amplitude of the second wave is the maximum while that of the first wave is the smallest. These features suggest that the prevailing pandemic wave is not periodic. Hence, the MWSIR model is further modified in order to explain its nonperiodic behavior. One can remember that transformation of SWSIR to MWSIR is accomplished by introducing two parameters γ and τ , where γ is the rate at which the recovered people become susceptible again, while τ is the time delay introduced during the estimation of susceptible population in terms of infected population.

In order to capture this specific nonperiodic feature in the most simple way, the value of τ , α and γ have been split up in the MWSIR model as shown below,

where τ 1 = 173, τ 2 = 440, τ 3 = 457, α 1 = 1.48 × 10 −7 , α 2 = 8.7 × 10 −7 , γ 1 = 0.055, and γ 2 = 0.08. The value of recovery rate, β = 1/14 per day, is kept fixed as it corresponds to the mean recovery period of COVID-19 infection or 1/β = 14 days. Values of those constants are estimated with respect to the single infected people, I (0) = 1. The numerical data obtained by solving the MWSIR model is drawn in red dashed line and compared with the daily new COVID-19 cases (seven-day moving average) in India drawn in blue sold line, which is shown in Fig. 6 . A very good agreement is found besides the fact that height of the first peak is higher and the width of the second peak is narrow in this estimation with respect to the actual values. The relation, α 2 > α 1 may attribute to the fact that transmission rate of new variants, like delta and omicron is higher than the primitive SARS-CoV-2 virus. In order to confirm the values of α 1 and α 2 geometrically both S(t) and I (t) are drawn in Fig. 7 . As the positions of peak of I (t) could be identified by the relation, This phenomenological study predicts that 5.5% of individuals infected during January 30, 2020 to April 14, 2021 become susceptible again for another infection. Similarly, 8.0% of individuals infected during later period become susceptible again for the next infection. These predictions are made in accordance to the estimated values of γ 1 = 0.055, and γ 2 = 0.08 in the different time domain. The difference, γ 1 = γ 2 , may account the effects of quarantine, isolation, vaccination and other preventive measures along with the presence of new variants like delta and omicron those are regarded as more infectious. However, as the solutions of the nonlinear equations are highly sensitive to the initial conditions, like S(0), I (0), and R(0), estimated values of the parameters, α, γ , and τ are likely to change accordingly with the change of initial values. Which on the other hand will change the predictions for obvious reason.

In order to study the behaviour of ongoing pandemic of COVID-19 in India mathematically, a modified version of SWSIR model termed as MWSIR model is pro-posed here which is found successful to reproduce the emergence of the multiple peaks in epidemiological wave prevailed in this country. The model is deigned in such a fashion that it is able to yield periodic as well as nonperiodic epidemiological waves with varying widths and heights by tuning its parameters. Equation 1a reveals that the sign ofṠ is determined by the factor α S(t) I (t) − γ I (t − τ ), when μ = 0. So it depends on the values of both S(t) and I (t) and as a result the opposite relations,Ṡ < 0 andṠ > 0 are being satisfied alternately with the increase of time. This property ultimately leads to oscillatory solutions.

As the SARS-CoV-2 virus responsible for the COVID-19 is highly infectious, the spreading of this disease cannot be controlled easily by means of simple preventive measures like quarantine, isolation, lockdown and even vaccination. In addition, emergence of super-spreading new variants, for example, delta and omicron makes the situation worse at later time. As a result, in due course, a sizable fraction of recovered people gets infected multiple times leading to multiple waves of the epidemiological infection. In order to explain this particular feature, the SIR model has been modified where a finite probability of the recovered people for becoming susceptible again after their recovery is taken into account. The waves of epidemiological evolution witnessed in India during January 30, 2020 to March 14, 2022 (775 days) have been reproduced satisfactorily in this formulation. SARS-CoV-2 virus was transmitted into India by a number of COVID-19 patients who entered this country mainly via air-travel mode through different airports of India continuously by the end of January 2020. Testing facility of COVID-19 infection was not adequate at that time. So the true number of infected individual for that time is not available. Therefore, this study is carried out with respect to the minimum number of infected individual, which means I (0) = 1. Estimated values of the parameters obtained in this study might be useful for imposing restrictions to control the pandemic in future. It will become helpful even for making accurate forecast and prediction.

However, in this simple model no time delay for the incubation period of the pathogens within the human body is considered, and the effect of quarantine, isolation, vaccination and other macroscopic measures are not take into account. So, it is expected that more accurate prediction can be made with the help of this model by accommodating the effect of quarantine, isolation, vaccination, lock-down and other factors in its modified version. Which means that the discrepancy between theoretical estimation and actual daily new cases can be removed further by formulating multi-wave SEIR, SIQR, SEIQR and other hybrid models. Obviously, this model can be employed for examining the characteristics of epidemiological waves found in other countries as well.

Funding The authors have not disclosed any funding.

The dataset that supports the findings of this study is available within the article.

Conflict of interest Authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the study reported in this article. Authors accept journal's ethical policies.

Contribution to mathematical theory of epidemics

The mathematics of infectious diseases

Mathematical epidemiology of infectious diseases

Modeling infectious diseases in humans and animals

An introduction to mathematical epidemiology

Epidemic modelling: an introduction

Effects of quarantine in six endemic models for infectious diseases

A SEIQR model for pandemic influenza and its parameter identification

Periodicity in epidemiological models

Applied mathematical ecology

Endemic models with arbitrarily distributed periods of infection, II: fast disease dynamics and permanent recovery

The construction of next-generation matrices for compartmental epidemic models

Estimation of COVID-19 outbreak size in Harbin

Using awareness to Z-control a SEIR model with overexposure: insights on Covid-19 pandemic

Analysis and forecast of COVID-19 spreading in China, Italy and France

Predicting the evolution of the COVID-19 epidemic with the A-SIR model: Lombardy, Italy and Sao Paulo State, Brazil

Understanding COVID-19 nonlinear multi-scale dynamic spreading in Italy

A Hybrid Compartmental Model with a Case Study of COVID-19 in

A stochastic SEIHR model for COVID-19 data fluctuations

COVID-19 pandemic in India: a mathematical model study

Dynamics of COVID-19 transmission with comorbidity: a data driven modelling based approach

Mathematical model of COVID-19 with comorbidity and controlling using non-pharmaceutical interventions and vaccination

Solving ordinary differential equations I