COVID-19: Development of a robust mathematical model and simulation package with consideration for ageing population and time delay for control action and resusceptibility.

The current global health emergency triggered by the pandemic COVID-19 is one of the greatest challenges we face in this generation. Computational simulations have played an important role to predict the development of the current pandemic. Such simulations enable early indications on the future projections of the pandemic and is useful to estimate the efficiency of control action in the battle against the SARS-CoV-2 virus. The SEIR model is a well-known method used in computational simulations of infectious viral diseases and it has been widely used to model other epidemics such as Ebola, SARS, MERS, and influenza A. This paper presents a modified SEIRS model with additional exit conditions in the form of death rates and resusceptibility, where we can tune the exit conditions in the model to extend prediction on the current projections of the pandemic into three possible outcomes; death, recovery, and recovery with a possibility of resusceptibility. The model also considers specific information such as ageing factor of the population, time delay on the development of the pandemic due to control action measures, as well as resusceptibility with temporal immune response. Owing to huge variations in clinical symptoms exhibited by COVID-19, the proposed model aims to reflect better on the current scenario and case data reported, such that the spread of the disease and the efficiency of the control action taken can be better understood. The model is verified using two case studies based on the real-world data in South Korea and Northern Ireland.

Introduction

The coronavirus disease COVID-19 is a respiratory infection disease caused by the novel coronavirus, SARS-CoV-2. The first COVID-19 outbreak was reported in Wuhan of Hubei Province, China at the end of December 2019. Within just two months, the disease has rapidly spread across the world and it has been declared a global pandemic in early March 2020. As of 20th April 2020, the virus has affected close to 2.5 million people with approximately 170,000 confirmed deaths across at least 184 countries [1].

The symptoms caused by the SARS-CoV-2 virus have large variations with most people only experiencing mild to moderate respiratory illnesses and only a smaller group of people would develop complications of respiratory failure or acute respiratory distress syndrome. Based on clinical data reported from Wuhan where the outbreak began, elderly patients have been identified to have higher odds to experience severe symptoms with higher mortality rate compared to people of younger age [2]. Study also shows that up to approximately 80% of the people infected with SARS-CoV-2 are asymptomatic carriers, i.e. they experience no or mild symptoms but are still able to transport the virus to others [3]. This has caused the detection and containment of the SARS-CoV-2 virus to be much more challenging. As a result, social distancing has been widely implemented in many countries worldwide to slow down the transportation of the virus through minimised human-to-human contact. Individuals who have recovered from COVID-19 after experiencing mild or moderate symptoms are more likely to develop temporary resistance towards the virus and are unlikely to experience severe respiratory illnesses [4]. However, in rare occasions, there have been clinical findings showing that patients who have recovered from the disease have been tested positive again. For instance, in February 2020, a patient in Osaka, Japan, has been tested positive towards the SARS-CoV-2 a few days after being discharged from the hospital for treatment with the disease [5]. Due to very limited knowledge on the immune response of humans to this novel virus, the possibility of reinfection cannot be ruled out at the moment.

Mathematical modelling and computational simulations have played important roles in describing the dynamics of infectious diseases using nonlinear systems so that their risks could be better understood [6], [7], [8], [9], [10], [11]. Most notably, the SEIR (Susceptible–Exposed–Infected–Removed/Recovered) model has been widely reported during the past decades in quantitative modelling studies of infectious epidemics/pandemics, such as the severe acute respiratory syndrome (SARS) in 2002 [12], the influenza A (pH1N1) pandemic in 2009 [13], the Middle East respiratory syndrome (MERS) pandemic in 2013 [14], as well as the latest Ebola outbreak in 2018 [15]. The SEIR model represents a typical infectious epidemic disease using four distinctive phases; susceptible (S) represents the population that has yet to be infected by the virus; exposed (E) represents the number of individuals exposed to the virus; infected (I) models the number of people infected, have demonstrated symptoms, and are able to spread the virus to the people in the S compartment; and lastly, recovered/removed (R) models the number of people who have recovered and are assumed to have immune response towards the virus [16], [17], [18], [19]. Thus, based on the model, the S compartment will slowly deplete as the outbreak prolongs further, and the virus will eventually die out due to insufficient population within the S compartment. This compartmental modelling method allows transport of population from one compartment to another, where the disease transmission rates with respect to time can be simulated.

In this work, we propose a modified SEIRS model based on the Kermack–McKendrick model [20] with consideration for time delay and resusceptibility to the virus after recovery. In this model, the probability of a recovered patient to be reinfected with SARS-CoV-2 is taken into consideration to predict the future projection of COVID-19 cases. Resusceptibility is one of the crucial keys that could possibility lead to a prolonged COVID-19 pandemic. Time delay in the control action representing the time taken for the authorities to act on the virus and also the duration of short-term immunity after recovery, which may lead to resusceptibility, are also considered in the model. The time delay factor is applied to enhance the robustness and accuracy of the model and simulation, hence to better reflect on the timely situation with specific measures to control the transmission of the disease. The consideration for resusceptibility with time delay is an important highlight in this paper as it has rarely been considered in the SEIR models reported thus far [18], [21], [22]. Other than that, we also included information such as demographic details for the ageing population, who seem to experience a higher fatality rate due to COVID-19 [23].

This paper is organised as follows: Section 2 introduces the mathematical model; Section 3 presents the theoretical proofs for the positivity, boundedness, and stability of the model, as well as describes the model with time delay factors for the control action and potential resusceptibility; Section 4 verifies the proposed model and also provides some extensive results and discussions on predictions using the model through case studies based on real-world data; and Section 5 concludes the paper. Appendix A presents the design of the simulation package using the MATLAB/Simulink environment.

Modelling COVID-19 using modified SEIRS

First, let us consider the modified SEIRS model system below, where , and represent the susceptible, exposed, infected, recovered/removed, and deaths compartments, respectively. It is established that , where is the total stock population. The constant is the birth rate in the overall population and is the death rate due to conditions other than COVID-19. The parameter is the rate of transmission per S-I contact, is the rate of which an exposed person becomes infected, and is the recovery rate. Therefore, the incubation and recovery times are and , respectively. The constant is the efficiency of the control action to reduce the infection rate and to flatten the curve. It has a direct effect on the basic reproduction number , which will be explained further later in this paper.

The parameter comes into effect in the worst case scenario where the patient does not recover from the virus. We model the fatality rate with influence of the fraction of elderly population (above 65 years of age) within the community , where the percentages of nonelderly and elderly who recovered are and , respectively. The time spent hospitalised or infected in fatal cases is . In this paper, we establish that , assuming that patients spend the same amount of time hospitalised or infected, whether they recover from the virus or not.

The function represents the resusceptible stock, which can be computed from the recovered population using where is the percentage of the recovered population who are resusceptible to the virus. The number of actual recovered/removed cases with permanent immunity can then be written using In an ideal situation, population who recovered develop permanent immunity against the virus, i.e. . As a result, (6) becomes and subsequently, (7) becomes .

Positivity, boundedness, and equilibrium analysis of the model

Positivity of the solutions

The solutions to all subpopulations in the system (1) – (5) are nonnegative for all time given any finite nonnegative initial conditions of .

Firstly, it is established that all subpopulations defined by the system (1)–(5) are continuously differentiable. As such, if all subpopulations have nonnegative initial conditions, and that if any of the subpopulations is zero at time , its derivative is nonnegative by inspection. Assume that , , and at time instant . Then, we can rewrite (1) using where we can establish that and hence, is nonnegative for all time . Next, assume that , , and at time instant . We can rewrite (2) using so that and hence, is nonnegative for all time . Assume that , , and at time instant . Eq. (3) then becomes where we can establish that and hence, is nonnegative for all time . Assume that , , and at time instant . Eq. (4) can be rewritten using so that and hence, is nonnegative for all time . Finally, assume that , , and at time instant . We can then rewrite (5) using where we can establish that and hence, is nonnegative for all time .

It can be seen that since none of the subpopulations would have a negative derivative at any time instant of when all other subpopulations are nonnegative, then it can be concluded that all subpopulations are nonnegative for all time . As a result, given that , the stock population is also nonnegative for all time . Hence, the proof is complete.

Boundedness of the solutions

The stock population is finitely upperbounded for any nonnegative initial conditions.

The dynamics of the stock population can be written using Assuming that , and since Lemma 1 has established that all subpopulations are nonnegative and given that all parameters are assumed to be positive, then (8) becomes We can then deduce that where an integration of the inequality (10) yields for all time . As a result, the stock population is finitely upperbounded and hence, the proof is complete.

Disease-free equilibrium

The disease-free equilibrium is locally asymptotically stable if the basic reproduction number .

The disease-free equilibrium can be obtained by equating equations (1)–(4) to zero, hence satisfying of which the disease-free equilibrium is given by . Eq. (5) can be removed from this analysis without loss of generality as other equations do not depend on it. It can then be shown that the Jacobian for (1)–(4) at is written using

The characteristic equation can subsequently be obtained by subtracting from the diagonal elements and then computing the determinant, which yields where The first two eigenvalues in (16) can be easily computed to be , and that they are negative. As for the remaining eigenvalues, they can be found by solving the quadratic in (17), which can be expended and represented using where

For the disease-free equilibrium to be stable, i.e. all eigenvalues are negative, it is required that

As such, the basic reproduction number with the control action is defined using where for a disease-free system that is locally asymptotically stable, we need to ensure that while an unstable would translate to . Hence, the proof is complete.

Should there be no control action taken, i.e. , then the basic reproduction number in (18) becomes which is similar to other models found in the literature [6], [9], [10]. See [24] and the references within for a brief study on using control theory to fight COVID-19.

Endemic equilibrium

The endemic equilibrium is locally asymptotically stable if the basic reproduction number .

Let us assume that the system (1)–(4), which other than the disease-free equilibrium has a unique equilibrium at such that It can be seen that the model has no positive endemic equilibrium if since , and would be negative, which indicate an unrealistic biological system. If , then we would have the disease-free equilibrium discussed earlier in Section 3.3. Hence, for a positive endemic equilibrium system, we would require that .

For the stability analysis of the endemic equilibrium, we use the Jacobian for (1)–(4) at , which is written using where . The characteristic equation can subsequently be obtained by subtracting from the diagonal elements and then computing the determinant, which yields where

Assuming that   , then the characteristic equation in (24) becomes It is obvious that the first eigenvalue is . It can also be shown that in order for the remaining eigenvalues to be negative such that the endemic equilibrium is locally asymptotically stable, i.e. , and , we would require that while an unstable would translate to . Hence, the proof is complete.

For the remaining of this paper where we verify the model and perform predictions using real-world data in Section 4, and also for the model used in the simulation package presented in Appendix A, we assume a closed population with negligible birth and death rates, i.e. . Time delay factors are also considered for the control action taken as well as resusceptibility. As a result, the system (1)–(5) becomes and hence, the basic reproduction number in (18) becomes

The time delay indicates the time taken for the control action to take effect in flattening the infection curve, where represents the time to initiate the control action after the first confirmed case at , and represents the time after the control action has been initiated but before the effects are evidenced in the outputs of the system. In practical scenarios, can be used to model the delay for the population to effectively respond to the rules introduced by the control action, such as social distancing, self-isolation, and lockdown.

The function represents the resusceptible stock with the consideration for temporal immunity, of which we can then rewrite (6) using where the time delay represents the duration of temporal immune response of the recovered population. Hence, we can also update the expression for the number of recovered cases introduced in (7) using In an ideal situation where population who recovered develop permanent immunity against the virus, and . As a result, (32) becomes and subsequently, (33) can be rewritten as .

The block diagram of the proposed SEIRS model with time delay is shown in Fig. 1. The system with time delay is assumed to be stable and will exhibit similar disease-free equilibrium and endemic equilibrium properties as the system without time delay provided that the time delay parameters are nonnegative, i.e. . Detailed discussions on the theoretical stability analysis of SEIR and similar epidemic models with time delay can be found in studies such as [19], [25], [26], [27].

Fig. 1

The block diagram of the proposed SEIRS system used in the simulation package.

Case Studies

Case Study 1: Verification using data in South Korea

South Korea is used as a case study due to the amount of data available given that it is one of the first few countries to be directly affected by the virus outside of China. Its first confirmed case was reported on 20th January 2020 [28]. The other reason is that South Korea is also one of the very few countries that managed to effectively flatten the infection curve and it has set itself apart from others in leading the fight against COVID-19, at least for the moment. For example, the country started vigorous testing among its population with contact tracing, especially those of confirmed and suspected cases during the early stage of the epidemic. The government accomplished this by maintaining a public database keeping track of mobile phones, credit cards, and other data of patients who tested positive [29]. Also, on 16th March 2020, the authorities in the country began screening every person, both domestic and international, who arrived at its airports.

As of 20th April 2020, there have been 10,674 confirmed cases and 236 fatalities in South Korea [1]. As a result, we used the following parameters for our simulation. First, we assumed that the population of South Korea to be approximately with an elderly population of about 15% () [30]. We then set the percentage of recovery to be 98% () for the general public [1] and a fatality rate of () for the elderly [31]. We then assumed the incubation time and recovery time to be , and , respectively in accordance with [32]. The basic reproduction number was set to (95% CI: 4.97–5.23) based on the early growth-rate of the epidemic in South Korea. The initial infected and exposed cases were assumed to be and , respectively. See Table 1. Fig. 2a shows the results of the initial fitting of the model based on the data in South Korea while Fig. 2b shows the projections of the model when no control action is taken. There are some minor discrepancies between the modelled values and the real-world data during the initial stage of the simulation as seen in Fig. 2a. This is absolutely reasonable and acceptable while modelling an actual epidemic as most countries are still coming to terms with the virus during the first month and the data do not usually represent the actual number of cases due to lack of testing for confirmed cases.

Table 1

Initial parameters used to fit the model to the data in South Korea.

Parameter	Value
Stock population, N	51.5×10^6
Fraction of elderly population, N_old	0.15
Percentage of recovery, κ	0.98
Fatality rate for elderly, 1−κ_old	0.08
Incubation time, τ_inc	5.1 days
Recovery time, τ_rec	18.8 days
Basic reproduction number, R_0	5.1 (95% CI: 4.97–5.23)
Initial infected cases, I(0)	4
Initial exposed cases, E(0)	80

Fig. 2

Subfigure (a) shows the initial fitting of the model onto the data in South Korea and subfigure (b) shows the projections of the model when no control action is taken.

Once we have the initial fitting of the model, we introduced control action in line with the mitigation and preventive measures taken by the government. Due to the aforementioned vigorous testing, contact tracing, and isolation efforts taken, we assumed that the control action has an efficiency of 88% (). As a result, the reproduction number could be reduced to 0.61. We also assumed that there was a time delay of 30 days since the first confirmed case before the control action was introduced () and a further delay of approximately 13 days before the control action could be properly executed in the community (). Fig. 3 shows the simulation results; Fig. 3a shows that the trajectory of the modelled infected and death cases matches the real-world data after the control action was introduced. Fig. 3b shows the simulation results until the model stabilises assuming no subsequent control action being taken to further reduce the reproduction number. Comparing Fig. 3b with Fig. 2b, the peak of the number of infected cases could be reduced from about 19 million cases to about 7500 cases.

Fig. 3

The projections of the model onto the data in South Korea when control action with an efficiency of 88% is taken. Subfigure (a) shows the projections during the first 100 days while subfigure (b) shows the projections until the system achieves stability.

Subsequently, beginning 20th March 2020, stronger infectious disease control measures for travellers coming from overseas were enforced, where all Koreans and foreigners with residence in Korea arriving from all countries would be subject to self-quarantine for 14 days upon entry. All short-term travellers will also be ordered to self-quarantine with exceptions for some limited special cases [33]. Around the same time, the Korea Centres for Disease Control and Prevention (KCDC) also started advising all people in the country to observe an “Enhanced Social Distancing Campaign” [34]. Inducing these control actions into the model produces the results shown in Fig. 4. The second control action adds another efficiency of 50%, hence bringing the reproduction number further down to .

Fig. 4

The projections of the model onto the data in South Korea when a second control action with an efficiency of 50% is taken. Subfigure (a) shows the projections during the first 100 days while subfigure (b) shows the projections up till 300 days.

Simulation with resusceptibility

One of the many uncertainties about COVID-19 is whether patients who have recovered from the virus will be reinfected in the future. There have been reports in the news that patients who recovered from the virus were tested positive for a second time after being cleared of the virus [5], [35], [36]. On the other hand, most health authorities believe that patients who recovered may develop an immunity towards the virus. However, it is not sure if the said immunity is short-term or long-term. Hence, further research is required to provide clinical proofs to this hypothesis.

As such, we repeated the simulation for the Case Study on South Korea without control action, but with the inclusion of a possibility of resusceptibility. Here, we assumed that 1% of the patients who recovered are resusceptible towards the virus (), where the patients develop temporal immunity of days, respectively after recovering from the initial infection of the virus. Fig. 5 shows the simulation results, where the first infection spikes shown in all subfigures are synonymous to the result presented in Fig. 2b. The subsequent infection spikes are the result of resusceptibility, depending on the period of temporal immune response. The results show new surges in infection cases after the specific in each case, which diminish over time as more people develop immunity towards the virus. Interestingly, for the result shown in Fig. 5d where days, it could also be used to reflect on the situation where the virus may exhibit similar characteristics as the seasonal flu or the pandemic influenza A (pH1N1) that it is most likely active during certain seasons of the year, e.g. autumn/winter for the seasonal flu and spring/summer for the pH1N1, in which case an annual vaccine administration is necessary [37], [38].

Fig. 5

Trajectories for the infected and fatalities in South Korea due to resusceptibility where it is assumed that 1% of the recovered cases are reinfected after a time span of temporal immunity of (a) 0 day, (b) 30 days, (c) 90 days, and (d) 360 days, respectively. However, these results only apply assuming that if there is no control action being taken to flatten the curve.

Case Study 2: Prediction using data in Northern Ireland

Given the location of which this research is based, data in Northern Ireland is used for prediction study of the model. The reports on confirmed and death cases are published daily since 24th March 2020 by the Northern Ireland Public Health Agency (PHA) via their Daily COVID-19 Surveillance Reports [39]. The first confirmed case was recorded on 27th February 2020, and as of 20th April 2020, the total number of confirmed cases stood at 2728 with 207 fatalities.

We used the parameters in Table 2 for the initial fitting of the model based on the data from PHA on the initial growth-rate of the epidemic in Northern Ireland. Fig. 6a shows the results of the initial fitting, with Fig. 6b depicting the projections of the infected and deaths if no control action is taken. We then simulated the model based on the control action carried out; most schools in Northern Ireland were closed beginning 18th March 2020 followed by a national lockdown initiated by the United Kingdom government on the 23rd March 2020. As such, we set (20 days) to correspond to said dates since the first confirmed case, and assuming that it took a further approximately 12 days for the public to respond to these control actions, i.e. , we obtained the simulation results as shown in Fig. 7. The results show that in order for the model to follow the projected trajectories of the data from PHA in Fig. 7a, the control action has to achieve an efficiency of about 45% (), which indicates that the reproduction number could be reduced to 2.75. Comparing Fig. 6b with Fig. 7b, the peak of the number of infected cases could be reduced from 650,000 cases to 350,000 cases.

Table 2

Initial parameters used to fit the model to the data in Northern Ireland.

Parameter	Value
Stock population, N	1.88×10^6
Fraction of elderly population, N_old	0.18
Percentage of recovery, κ	0.94
Fatality rate for elderly, 1−κ_old	0.12
Incubation time, τ_inc	5.1 days
Recovery time, τ_rec	18.8 days
Basic reproduction number, R_0	5.0 (95% CI: 4.85–5.15)
Initial infected cases, I(0)	3
Initial exposed cases, E(0)	60

Fig. 6

Subfigure (a) shows the initial fitting of the model onto the data in Northern Ireland and subfigure (b) shows the projections of the model when no control action is taken.

Fig. 7

The projections of the model onto the data in Northern Ireland when control action with an efficiency of 45% is taken. Subfigure (a) shows the projections during the first 70 days while subfigure (b) shows the projections until the system achieves stability.

Further control action to meet critical Care Capacity

However, based on the results shown in Fig. 7, it is essential to further flatten this curve due to the limit of about 330 critical care beds available in Northern Ireland (100 initial setup + 230 introduced by the newly built Nightingale hospital) [40]. According to the Intensive Care National Audit & Research Centre (ICNARC) with its “Report on 2249 patients critically ill with COVID-19” dated 4th April 2020, about 6% of those tested positive for the SARS-CoV-2 required critical care [41]. Meanwhile in Italy, as of 29th March 2020, up to 12% of all positive cases were admitted to the intensive care unit (ICU) [42]. As such, should we assume that approximately 10% of those tested positive in Northern Ireland would require ICU admission, then the peak of the infection curve should not exceed 3300 cases, i.e. more control actions have to be taken to reduce the peak of 350,000 cases as seen in Fig. 7b.

Therefore, on day 38, which is about one week after the infection curve started to flatten due to the first control action, a second control action was introduced into the model. This second control action also reflects on the announcement made by the United Kingdom government in early April 2020 to allow police officers to enforce social distancing measures. Assuming that this second control action results in a further efficiency of 66%, the reproduction number could then be reduced to 0.93, and that it would take another 7 days for the public to fully respond to the control action, we could observe the results as shown in Fig. 8. With the initiation of the second control action, it can be seen in Fig. 8a that the peak in the infection curve is now reduced to 3500 cases. As such, the critical care capacity should be able to meet the demand for treatment based on the same assumption made earlier in this section, where it is estimated about 10% of the infected cases are admitted to the ICU. Another observation worth noting is that the number of deaths in the worst case scenario has also been reduced to about 3000 cases. See Fig. 8b. This value echos the projection made by the government in Northern Ireland that COVID-19 could lead to 3000 deaths during the first wave [43].

Fig. 8

Subfigure (a) shows the initial projections of the infected and cumulative deaths curves in Northern Ireland after the second control action with an efficiency of 66% was introduced on day 38 of the model. Subfigure (b) shows the remaining projections of the model.

Conclusion

This paper has presented a robust model for COVID-19 based on a modified SEIRS method to include considerations for the ageing population, and time delay for control action as well as resusceptibility of the recovered population due to temporal immunity. Two case studies using real-world data were presented in this research; the first case for verification of the model based on the data in South Korea including a study on the possibility of resusceptibility of recovered population; and the second case for prediction study of the model using data and up-to-date control action and related events in Northern Ireland. The simulation results from the case studies have clearly shown that the time of which the control action is taken and also the time for the public to properly respond to such intervention measures are critical in helping to flatten the curve. Also, until a time where a vaccine is developed and made available to the general public, the possibility of resusceptibility, no matter how small, will lead to subsequent waves of infections in the future depending on the time of temporal immunity. A simulation package was developed using the MATLAB/Simulink environment to ease understanding on the spread of the virus as well as the efficiency that needs to be achieved by the control action in order to successfully flatten the infection curve to not overload the healthcare capacity.

Interesting future research and expansion of the model include but not limited to the predictions for the occupancy of ICU beds, the effects of easing control action on and hence, the time of which control action has to be reinstigated, as well as specific subregions analysis such as demographic information to model the transmission of the virus among subregions in the country to cater for population movements and travels.

CRediT authorship contribution statement

Kok Yew Ng: Methodology, Software, Formal analysis, Visualization, Writing - original draft. Meei Mei Gui: Conceptualization, Data curation, Formal analysis, Writing - review & editing.