Numerical modelling of coronavirus pandemic in Peru

Mathematical model

Mathematical models are those that show the different scenarios that their solutions present according to the conditions established in the parameters and/or initial conditions, to which they become dependent for determining the stability or instability of dynamic systems.

Dynamic system models are always approximated based on several modelling assumptions, and epidemiological models are no exceptions to this basic principle. The right level of approximation crucially depends on the purpose of the model; depending on that, some assumptions may be perfectly acceptable or plain wrong.

The equations that govern the dynamics of an epidemic form a system of nonlinear first-order ordinary differential equations, the solution of this system is obtained analytically (under certain conditions) or by numerical methods, such as the finite difference method, Euler method and the Runge Kutta method (Cellier and Kofman 2006; Nakamura 2002).

In general, a numerical epidemic model is not exact and has uncertainties. Furthermore, the dynamics of the epidemic does not follow the mathematical model, but it does depend on the behaviour of the population. However, numerical models are useful to simulate epidemic scenarios that help in decision-making, at the central, regional or local government level.

The simplest mathematical model that describes the behaviour of an epidemic is the SIR model (susceptible – infected – removed), which considers groups of populations from the different phases of the epidemic process. In this model, the removed population includes the recovered population and those who died due to the disease. The SIR model assumes a relatively small time scale, such that the number of births and deaths due to causes other than the epidemic is insignificant. The governing equations of the SIR model are (Bacaer et al. 2021; Brauer and Castillo-Chavez 2012; Martcheva 2015):

Where S(t) is the susceptible population, I(t) is the infected population, R(t) is the removed (recovered + deceased) population, n is the total population, β is the infection rate and γ is the recovery rate.

The susceptible population will become infected at a rate of β, therefore the population change will be −βSI/n. The infected group feeds on the removal of the susceptible group at a ratio βSI/n but decreases due to the recovered (or deceased) at a ratio −γI. The removed population increases due to the infected that recover (or die) at a rate γI.

The spread of an epidemic is quantified by the basic reproduction number R ₀. This dimensionless quantity is used to represent the average number of infected produced by a member of the infected population during the contagion period (Brauer and Castillo-Chavez 2012). For the SIR model, it is defined as:

If R ₀ > 1, the epidemic will spread, while if R ₀ < 1, the epidemic will eventually die out. In general, in the case of coronavirus, this parameter varies between 1.5 and 3.5.

Kermack and McKendrick (1927) investigated the SIR numerical model for epidemics, they developed the mathematical treatment of the problem. In our country (Perú), Lopez, Vidal, and Valdez (2015) have reviewed the variants of the SIR model in a didactic way.

A variation of the SIR numerical model arises when the deceased population is separated from the removed population, considering that the rate of change of the deceased population is proportional to the infected people. This new numerical model is called susceptible – infected – recovered – died (SIRD) (Osemwinyen and Diakhaby 2015).

When the population in quarantine Q(t), social isolation A(t) and the effect of vaccination ψS/n are incorporated into the numerical model, then a new numerical model is obtained: susceptible – isolated – infected – recovered – deceased – quarantined (SAIRDQ).

We are not taking into account a group of vaccinated people V(t) because they are within the recovered population R(t). The effect of vaccination is represented by the term ψS/n in Eqs. (5) and (8). In Figure 1, the effect of vaccination is represented by the loop from the compartment S(t) to R(t).

Figure 1:

Flow diagram of the dynamic evolution model of the epidemic: SAIRDQ with vaccination. The model parameters are explained in the text.

Covid-19 is well-known to have a significant incubation period (around 5–7 days); this parameter is essential to forecast the short-term evolution of the epidemic (a few weeks). This issue is not included in the model because we are concerned about the long-term evolution of the epidemic.

We are not taking into account the case of the reinfected population due to the new variants of Coronavirus (delta variant, mu variant, etc.). It has been reported that few recovered people have been reinfected after some months.

Covid-19 is well-known for a significant fraction of asymptomatic infected subjects. The ratio of asymptomatic to infected people has been the subject of speculation since the beginning of the pandemic and is probably 40% or more. In this compartmental model, we are taking into account the undiagnosed infective people (asymptomatic) within the infected population I(t). This makes the model correct from a conceptual point of view.

This big hidden reservoir is also important for the spread of the epidemics, asymptomatic carriers are known to contribute significantly to the circulation of the disease, they can cause other people to contract symptomatic Covid and possibly to die. However, the numerical modelling of two groups of infected populations (symptomatic and asymptomatic) would go beyond the scope of this article and would be the subject of new research.

Figure 1 shows the flow diagram of the SAIRDQ model. The non-linear differential equations that govern this model are:

Where S(t) is the susceptible population, A(t) is the isolated population taken from the susceptible group, I(t) is the infected population, R(t) is the recovered population, D(t) is the deceased population, Q(t) is the infected population in quarantine, n is the total population, α is the isolation dropout rate, β is the rate of infection, γ and γ ₂ are the rates of recovery, η and η ₂ are the isolation rates for isolation and μ and μ ₂ are the mortality rates from the disease. The parameter ψ is the vaccination rate. If the parameters α and η are cancelled, then the isolation and quarantine effects would be cancelled and the SIRD model would be obtained.

If the population behaviour were homogeneous with no changes, the epidemiological parameters would be constant and the corresponding curves would be very simple; this will generate only one epidemic wave. On the contrary, the population behaviour in Perú was variable and not homogeneous; therefore, the epidemiological parameters are not constant during the evolution of the pandemic and the corresponding curves were complex with the presence of a second epidemic wave.

The basic reproduction number R ₀ varies between 2.0 and 3.0 in the case of Peru. Munayco et al. (2020) found a value for the initial phase of the pandemic in Lima of R ₀ = 2.3, which can be extrapolated for Peru (Munayco et al. 2020). For the SAIRDQ model, the parameter R ₀ is defined as:

From Eq. (11), the basic reproduction number was estimated at R ₀ = 2.30.

Numerical model

Eqs. (5)–(10) form a nonlinear system of first-order ordinary differential equations. It is possible to conduct a qualitative analysis of this kind of ordinary differential equation system and to find an analytical solution to these equations under certain conditions. Many authors have investigated in this regard (Giordano et al. 2020; Kermack and McKendrick 1927). However, from a numerical or computational point of view, it is easier to obtain a numerical solution by applying the finite difference method, Euler method or Runge Kutta method (Carcione et al. 2020; Cellier and Kofman 2006).

The Euler method is only a first-order accurate method. It is hardly a surprise that the solution using a fixed time-step algorithm be very sensitive to this time step Δt, particularly given the unstable nature of some parts of the paths, where the linearized system has positive eigenvalues.

The ordinary differential Eqs. (5)–(10) can be discretized using the finite difference method:

Where k is the discrete-time variable and Δt is the computational time step. For this case, a value (Δt = 0.01 day) much less than the sampling interval of the time series (T _s = 1 day) has been chosen.

In this case, the results of the numerical model are very sensitive to the choice of the computational time step, the increase in the value of Δt causes a decrease in the values of I(t), R(t) and D(t). From the definition of the derivative, if the time step Δt → 0, then the number of iterations Nk → ∞. A good choice was Δt = 0.01 day.

The regime values (when derivative is zero) can not change with the integration step size (Δt). Therefore, S = S ₀ exp(R ₀(S − 1)), I = 0, Q = 0 and A = (η/α)S, independently on Δt.

A routine has been implemented in the Matlab programming language that performs the recursive process using the finite difference method. Previously, the initial conditions, the model constants (α, β, γ, η and μ), the time step Δt and the number of iterations Nk must be defined.

To start the recursive process, it is necessary to have the initial conditions of the system: S(1), A(1), I(1), R(1), D(1) and Q(1). The initial conditions have been chosen based on the reality of Peru. The total population is n = 32.5 millions, it is considered that at the beginning of the epidemic (March 6, 2020) there was only one infected people I(1) = 1, the initial susceptible population S(1) = N − I(1) = 32,499,999. The initial population in isolation, quarantined, recovered and deceased is null: A(1) = 0, R(1) = 0, D(1) = 0 and Q(1) = 0.

Inversion: iterative approximation method

The inversion method is used in geophysics to obtain the parameters that characterize a system from the observations and data. In seismology, the inversion method is used to obtain the seismic source parameters from the seismic recordings, tsunami waveforms and geodetic data (Jiménez et al. 2014, 2020). In this research, we have used the inversion method (iterative approximation method) to obtain the epidemiological parameters from the reported data.

The challenging point of this research is the choice of the values of the epidemiological parameters, which are variable or constant piece-wise (Table 1). This is the work of the epidemiologist specialist. However, the choice of these parameters is possible under certain empirical criteria.

Because there are many parameters to be calculated and the system would be unstable, in the beginning, we have fixed some parameters. For example the recovery rate γ = 1/7 from the data published in the literature, the initial vaccination rate ψ = 0, because the vaccination process began in 2021, etc.

Some parameters (β and η) are not constant (or they are constant piecewise), these parameters vary in time according to many factors, such as the population behaviour concerning the pandemic, sanitary measures given by the government, etc. The available data were divided into several separate sections, with boundaries at 90, 120, 196, 200, 305, etc., from the first day. These intervals are related to specific non-pharmaceutical interventions taken by the government of Peru, as explained in Table 2.

To estimate the epidemiological parameters, we used the data inversion method, in which the observed values (number of deceased people) are compared with the simulated ones and the residual is reduced until a threshold error is reached (Marinov et al. 2014). The alternative is to use the iterative approximation method, in which some system parameters are fixed and the rest are varied until a good correlation is obtained between the observed and simulated data (Figure 5, Table 1). A metric to evaluate the correlation would be the normalized variance, which is defined as:

Where obs(k) represents the observed variable (number de reported deceased people), sim(k) is the simulated variable (number of simulated deceased people), and Nk is the number of data. Figure 6 shows the sensitivity test to find the parameters: initial β and γ. Values that minimize the normalized variance are taken into account. We have obtained a normalized variance equal to 9.8 × 10⁻⁵. Therefore, the observed and simulated death curves almost overlap (Figure 5).

Figure 5:

According to the results of the SAIRDQ model, there will be around 200 thousand deaths at the end of the pandemic. The blue curve is the number of simulated infected per day and the red curve is the simulated cumulative number of deaths. The yellow line represent the deceased data. The day 1 is 06 March 2020.

Figure 6:

To the left: sensitivity test to calculate the value of β = 0.48305 day⁻¹. To the right: sensitivity test to calculate the value of γ = 0.14102 day⁻¹. The normalized variance is 9.8 × 10⁻⁵.

It is important to notice that we are using the number of reported deceased people rather than the number of reported infected subjects in the inversion process. Unfortunately, the number of reported infected people is always underestimated compared to the actual number of infected people. This happens because many subjects that do not show symptoms are not tested, therefore they escape the official report. Therefore, the inversion process is based uniquely on the deceased data.

The numerical model was run on an Intel Core i7 computer, using a Matlab-encoded program. The duration of the program execution took 3 s on average, for a time step Δt = 0.01 day, with a total of Nk = 60,000 iterations.

Effect of social isolation and quarantine

We have simulated the scenario with the relaxation of social isolation and quarantine from the end of the first wave of the epidemic. The death toll would be around 1.5 million people. In this case, the second wave would be much larger than the first one.

Figure 7(a) shows the evolution of the pandemic taking into account the finish of the quarantine on 28 February 2021 and then the conditions are relaxed (η = 0, η ₂ = 0). We noticed that the second wave would be a “tsunami”. The simulated death toll would be around 280 thousand.

Figure 7:

To the left: the quarantine is finished on 28 Feb 2021 and then the conditions are relaxed. The second wave would be a true “tsunami”, with more than 280 thousand deaths in the next 12 months. To the right: no vaccination was applied, the death toll would be more than 400 thousand people. Day one corresponds to March 6, 2020.

On the other hand, in the ideal case of total social isolation and quarantine (at 100%), the peak of infections would be 8 and only 4 people would lose their lives, that is, those infected would be only in the family environment of the first infected.

Due to the end of the national quarantine on July 1, 2020, and the relaxation conditions of the population behaviour concerning quarantine since before the finish of it, the parameter β and η were considered variable by sections, according to Table 1. Therefore, the conditions of the SAIRDQ numerical model have changed. In Figure 5, the base of the infection curve (in blue) has widened, which means we have a second epidemic wave.

One effect of the quarantine and social isolation is the widening and flattening of the epidemic curves and the existence of successive waves of the epidemic. This is good for providing valuable resilience time to the overloaded healthcare system.

Effect of the vaccination process

We have simulated the effect of no vaccination and the same conditions of parameters in Table 1. The results show a death toll of more than 400 thousand people and the presence of at least two waves of the epidemic. In this scenario would be more than two waves and the death toll would increase with time.

In the real case, the vaccination began on 8 February 2021 (day 340) with the highest age group and so on descending. In the first weeks, the vaccination rate was very slow, ψ = 12,000 on average, due to global vaccine shortages. From the day = 416 (25 April 2021), the vaccination rate has increased to ψ = 90,000 on average. Now, the vaccination rate is: ψ = 100,000 on average.

To date (December 2021), more than 21.5 million people have been vaccinated with the two doses of the vaccine. This represents 66% of the total Peruvian population (https://data.larepublica.pe/avance-vacunacion-covid-19-peru). It is possible that by the end of the year 2021, the entire population over 12 years old in Peru will be vaccinated.

Vaccination is the most important way to defeat the coronavirus epidemic. As a recommendation, we must continue in social isolation and family quarantine until the arrival of vaccines and vaccination of all the family members.

Estimation of percentage of undiagnosed infective people

We have conducted the inversion process taking into account only the number of deceased people, because this parameter is more confident than the number of infected people, due to the number of undiagnosed infective people (asymptomatic).

The number of reported infected subjects is always underestimated compared to the actual number of infected subjects; therefore, the difference between the simulated infected people and the reported infected people would be the number of undiagnosed people (Figure 8B). This big hidden reservoir is also very important for the spread of the coronavirus. An estimation of the percentage of asymptomatic people is (Eq. (19)):

Where: I _r is the number of reported infected people (2.2 million) and I _s is the number of simulated infected people (3.99 million). Therefore, according to Eq. (19), the asymptomatic group represents the 45% approximately. This result agrees with that of Ma et al. (2021), with a global percentage of asymptomatics of 40.5%.

Figure 8:

Comparison of epidemiological curves.

(A) Comparison of the reported and simulated deceased data which were tuned in order to minimize the normalized variance. (B) Comparison of reported and simulated infected people, the difference represents an estimation of the number of undiagnosed infective people.

The ratio of reported deceased people to reported infected subjects in Peru (around 10%) is much higher than in most other countries (typically 1–2%), which may be due, at least in part, to an underestimation of the actual number of infected people. Statistics miss the asymptomatic subjects.

About the new variants of the coronavirus

In the last months, new variants of the coronavirus have appeared. The emergence of new, more contagious variants, such as the Delta and Omicron variants, as well as the increasing evidence that immunity (either from vaccination or from illness) can fade over time, would probably make the model proposed in this research unsuitable to handle the next years of the pandemic. However, these new factors would not be important during the period considered in this research (2020–2021).

To take into account the new variants of coronavirus, the numerical model would change with the addition of a path from the R compartment back into the S compartment; however, this would be the subject of another research paper.