Inverse problem for adaptive SIR model: Application to COVID-19 in Latin America.

This work presents a method for solving an Adaptive Susceptible-Infected-Removed (A-SIR) epidemic model with time-dependent transmission and removal rates. Available COVID-19 data as of March 2021 are used for identifying the rates from an inverse problem. The estimated rates are used to solve the adaptive SIR system for the spread of the infectious disease. This method simultaneously solves the problem for the time-dependent rates and the unknown functions of the A-SIR system. Presented results show the spread of COVID-19 in the World, Argentina, Brazil, Colombia, Dominican Republic, and Honduras. Comparisons of the reported affected by the disease individuals from the available real data and the values obtained with the A-SIR model demonstrate how well the model simulates the dynamic of the infectious disease.

Introduction

The spread of epidemic diseases has attracted attention for centuries, due to the importance of predicting their dynamics (Allman & Rhodes, 2003; Anderson & May 1991; Bauch et al., 2003; Bellomo & Preziosi, 1995; Diekmann & Heesterbeek, 2000). For example, Daniel Bernoulli analyzed smallpox morbidity in 1766, according to (Hethcote, 2000). Common infectious diseases include measles, malaria, varicella, HIV, Ebola, and SARS. There is no vaccination for some of the infectious diseases; however, there exist preventive practices for many of them, see (Kabir, Kuga, & Tanimo, 2019). The purpose of the developed mathematical models is to predict the future of the epidemic spread. Simulation results assist governments in making decisions on how to deal with the disease (Jin & Jia, 2020). It must be noted that mathematical models have limitations because they work under certain assumptions. For example, the reported COVID-19 epidemic data does not include all infectious cases. Decisions on governmental restrictions and vaccinations cannot sometimes be known and predicted.

The SIR (Susceptibles–Infected–Recovered) model origins can be traced back to the work of Kermack and McKendrick (Kermack & McKendrick, 1927) in 1921. This theory attempts to explain the rapid rise and fall in the number of infected people with a contagious illness in a closed population over time. This model is the basis of all current modeling of the dynamics and evolution of infectious diseases (Azam et al., 2020; Murray, 1993; Smith & Moore, 2004; Takeuchi, Iwasa, & Sato, 2007; Tanimoto, 2015, 2018; Tarantola, 2005).

There exist other models considering additional events and features, such as latent period, vaccine effect, reinfection, non-constant population, etc. For example, the SEIR (E stands for exposed) model is successfully used for COVID-19 dynamics prediction in (Kabir & Tanimo, 2020). Other models include SIUR (U stands for unreported cases), but not limited to, used in Liu and collaborators (Liu, Magal, Seydi, & Webb, 2020, Liu, Magal, & Webb, 2021). This model is used for predicting the number of reported and unreported cases for the COVID-19 epidemics in several countries. There is a number of recent publications on COVID-19 aiming to provide insight and understanding the trends of the disease (Ajbar, Alqahtani, & Boumaza, 2021; Griette, Magal, & Seydi, 2020; Griette & Magal, 2021; Kucharski, Russell, & Diamond, 2020; Li et al., 2020; Lin et al., 2020; Lobo et al., 2020; Nishiura, Linton, & Akhmetzhanov, 2020; Pereira, Schimit, & Bezerra, 2021; Roosa et al., 2020; Shereen, Khan, Kazmi, Bashir, & Siddique, 2020; Wacker & Schlüter, 2020).

The SIR model is a highly nonlinear dynamical system, particularly so in the case of COVID-19, due to its diverse characteristics. The rates of transmission and removal in the case of the COVID-19 depend on the evolution of the epidemic disease over time, see (Acedo, Morano, Santonja, & Villanueva, 2016; Liu, Wei, & Zhang, 2019). Every new epidemics must be studied and mathematical models constructed to answer important questions on handling the disease (Coronel, Huancas, & Sepúlveda, 2019). A method for the identification of the time-dependent transmission rate is proposed and used by other researchers, see the recent work (Demongeot, Griette, & Magal, 2020).

Recovering coefficients and other parameters in differential equations from overposed data is an example of inverse problem. In general, a problem is inverse if the values of some model parameter(s) must be obtained from the observed data. The definition of inverse problem involving differential equation(s), according to (Bellomo & Preziosi, 1995), is: “An initial-boundary-value problem is inverse if some information on the initial and/or boundary conditions needed for solution or/and on the parameters that characterize the model are missing and are replaced by suitable information on the solution of the mathematical problem.”

Part of the inverse problems are so-called parameter identification problems: adjusting the parameters to reproduce measured data. Since the number of infected and recovered persons always include random errors, a method for smoothing the data in order to evade the instability provoked by the pollution of the data was introduced in (Marinov, Marinova, Omojola, & Jackson, 2014) and later modified in (Marinov & Marinova, 2020). The idea of the method is to replace the incorrect problem with a well-posed problem for minimization of quadratic functional of the original equations. This way of identifying the unknown coefficients employs the values of infected and recovered persons, avoiding any statistical methods or artificial assumption about parameters of the model.

The present work utilizes a novel inverse problem approach to the time-dependent transmission and removal rates identification in the A-SIR (Adaptive SIR) model. The spread of the epidemics can be predicted by using the estimated parameter values to solve the A-SIR system and obtain the dynamics of the infectious disease spread. If conditions change, then the predictions may no longer be accurate; hence, adjustments will be required. We apply the method for the A-SIR model to investigate the COVID-19 epidemic for countries in Latin America. This gives insight into how well the method simulates the spread.

The A-SIR model for the spread of an infectious disease

The SIR model for the spread of an infectious disease was proposed by Kermack and McKendrick (Kermack & McKendrick, 1927). SIR stands for Susceptible, Infectious and Recovered. Numbers in each group are functions of time denoted by: S(t) – susceptible who can catch the disease, I(t) – infectives who have the disease and can transmit it, and R(t) – removed individuals who have either had the disease, or are recovered, immune or isolated until recovered. The SIR model assumes that the removed individuals are no longer susceptible nor infectious. The total population N = S(t) + I(t) + R(t) is considered constant in the classical model.

Specific details on the SIR model can be found in (Murray, 1993) where the classical model assumes that the rates are constants. The rate of change of S(t) is − βS(t)I(t), where β > 0. The infectious individuals leave I(t) with rate γ and they move directly into the R(t) group. The number of cases for recovered from COVID-19 individuals who are re-infected at the present moment is very limited and the rate cannot be estimated; thus, this possibility in not taken into account.

However, in the case of a pandemic, the rates may vary in time. Multiple factors may cause the rates to change over time. In the case of COVID-19, examples include social distancing, restrictions imposed by governments, and preventive treatments. In this work, the model considers time-dependent rates, namely β = β(t) and γ = γ(t). The adaptive SIR model for the Susceptible-Infected-Removed (Marinov & Marinova, 2020; dos Santos, Almeida, & de Moura, 2021), called here A-SIR model. A-SIR is a variation of the SIR model, consisting of the following three equations:

Fig. 1 shows the diagram for the A-SIR model, corresponding to the system of equations (1), (2), (3).

Fig. 1

The A-SIR epidemic model.

A fundamental question is under what circumstances an epidemic occurs. This happens when the rate of change of the class I(t) is positive, i.e., the right-hand side of equation (2) satisfies β(t)S(t)I(t) − γ(t)I(t) > 0, which is equivalent to

The parameter R(t) is called effective reproduction number (or effective reproduction rate or ratio) for a given disease. An epidemic occurs if an infective individual introduced into a population of susceptible individuals S(t) infects on average more than one other person. When the fraction of the population that is immune increases (because of vaccination or because of recovering from the disease) so that R(t) < 1, then herd immunity takes place and the number of new cases in the population decreases to zero.

Another important characteristic is the basic reproduction number (ratio, rate) R0(t), which is the ratio of the rates β(t) and γ(t) over time multiplied by the size of the total population N:

Normally, the basic reproduction number R0 is a constant. However, in the case of COVID-19, it depends on other factors, such as restrictions introduced by governments and societal behaviour. For this reason, we consider the basic reproduction number as a function of time.

The initial-value problem consisting of the system (1)–(3), with coefficients β(t) and γ(t) known, along with initial conditions derived from the given data, constitutes the direct problem. The inverse problem consists of determining the coefficients from the available data.

The inverse problem for the time-dependent rates and unknown functions

Let us assume that the values of S(t) and I(t) are known from the available data at two time moments, namely an initial time moment t = T and a final time moment t = T,

The original SIR model assumes that the transmission and removal rates are constants. Equations (1), (2), and the initial conditions (6) allow the determination of I(t) and S(t), if the coefficients β and γ are known constants. In this case, the terminal conditions (7) may not be satisfied exactly, because the problem becomes overdetermined. A method for solving the inverse problem for the classical SIR model is presented in (Marinov et al., 2014), where it is assumed that the problem is correctly posed according to Tikhonov (Tikhonov & Arsenin, 1974).

The inverse problem approach to solving the adaptive SIR model (A-SIR) considers the time-dependent transmission and removal rates as piece-wise constant (step) functions of time

A dataset of values at given time moments t1, t2, …is assumed to be available: D = {(S(t), I(t)), k = 1, 2, …}.

Then, the constants β and γ can be estimated by solving the inverse sub-problems (1), (2), (6), (7) given that S = S(t), S = S(t), I = I(t), and I = I(t), for k = 2, 3, …

If β and γ are constant and unknown, the general solution of the system (1), (2) depends on four constants – two constants from the integration and two constants from the unknown coefficients. There are four equations in (6), (7); hence, the problem for identifying simultaneously the coefficients β, γ, and the functions S(t) and I(t) is well-posed under the assumption that a solution exists. A preliminary version of the algorithm for solving this inverse problem using early COVID-19 data is given in (Marinov & Marinova, 2020).

The described algorithm is applicable to a dataset D that includes exact values of all infectious individuals, which is far from being the case in the COVID-19 epidemic. Available data contains random noise and errors, which cause strong oscillations in the obtained solution of the formulated above inverse problem, seen also in the results of our simulations. Therefore, an algorithm must address these inaccuracies. Next section is devoted to presenting a method for solving the inverse problem for time-dependent rates and unknown functions, which handles data inaccuracies.

Solving the inverse problem

The solution method described here permits the recovery of time-varying coefficients, which takes into account the non-autonomous nature of an epidemic. In other words, the population adapts to the presence of the disease and the global behaviour of the epidemic, which is changing over time.

Minimization problem using data values at two time moments

First, we start with the sub-problem using available data values at two time moments. Since we are concerned with the numerical solution of the system (1)–(3), we seek approximations of the functions S(t), I(t), and R(t) at the discrete set of points {t0, t1, …, t} in the interval [T, T], where T is the initial time moment, T is the final time moment, and n > 1, see Fig. 2.

Fig. 2

The grid of equidistant points t = T + kτ, k = 0, 1, …, n.

The time step is defined to be , and the equidistant nodes are given by t = T + kτ, k = 0, 1, …, n.

Since the inverse problem is non-linear, estimating the values S and I requires iterations. Let and denote values from a previous iteration, and , for k = 1, 2, …, n.

Now consider the second order discretization of equations (1), (2) on the regular grid corresponding to Fig. 2:

Note that equation (9) also involve the unknown rates β and γ on the right-hand side, which can be determined from the available data. Therefore, the problem is inverse and requires a special treatment.

Next, we describe our approach to solving the inverse problem. Let ε and δ be the residuals of equation (9), namely

A solution of the inverse problem is given by the set of values minimizing function

Assuming approximate values of β and γ are known, the necessary conditions for minimization of the function Φ with respect to its arguments S and I yield the following linear difference equationsfor k = 1, …, n − 1. After adding initial and terminal conditions (6) and (7) to equations (11), (12), the obtained linear system of 2(n + 1) equations for the unknowns (S0, S1, …, S) and (I0, I1, …, I) becomes well-posed. This way of linearization allows using a matrix with constant elements on the left-hand side of the system of equations; thus, producing an efficient iterative process of approximating the solution.

Now assuming S and I are known, we derive explicit formulas for the rates β and γ from the necessary conditions for minimization of the function Φ in equation (10) with respect to β and γ:where

Details on the derivation of equations (11), (12), (13) are given in (Marinov & Marinova, 2020).

Minimization problem using the entire dataset

Let the number of infectious individuals at some time moments ν1, ν2, …, ν be known and given byfor l = 1, 2, …, m. The values of the coefficients β and γ are considered unknown.

Assume that for every 1 ≤ l ≤ m, there exists an index k, such that , i.e., the set of time moments {ν1, ν2, …, ν}, is a subset of the set of nodes {t0, t1, …, t}. Let χ denote a value σ such that: if there exists k ∈ {1, …, n} with t = ν, 1 ≤ l ≤ m, then χ = σ and μ > 0; otherwise χ = 0 and μ = 0, where μ is the weight of the residual of equation (14) in the function Ψ.

Similarly to Subsection 4.1, a solution is obtained as a result of minimization of a function. The function is chosen to bewhere ε and δ are the residuals of equation (9).

The necessary conditions for minimizations of the function Ψ with respect to S and I (k = 1, 2, …, n − 1) give the equations

After adding the initial and terminal conditions (6), (7), the linear system (16), (17) for the unknowns (S0, S1, …, S) and (I0, I1, …, I) becomes well-posed. The equations for β and γ in this case are identical to equation (13).

Time-dependent parameters β, γ, and R

Equation (14) is the dataset with the number of infected individuals for every day for a period of m days. The set is divided into subsets of fixed length of P days, as shown in Fig. 3. In order to approximate the time-dependent transmission and removal rates β and γ, Algorithm 1 is applied to every sub-interval [k − P + 1, k], k = P, P + 1, …, m. After computing the values of β and γ, the effective reproduction rate R = βS/γ is evaluated on every sub-interval.

Fig. 3

The subsets of fixed length of P days for identifying β(t) and γ(t).

Algorithm 1 shows the iterative procedure for solving the system (16), (17), (6), (7) along with obtaining the values of the transmission and removal rates from equation (13).

Iterative algorithm for solving the inverse problem.

A-SIR simulations using COVID-19 data over nine month period

Obtained numerical results from the simulations using available reported data for the COVID-19 pandemic over nine month period are reported here. We consider the World and selected countries in Latin America.

For visualizing purposes, several parameters are monitored and analyzed for better understanding of the infectious disease dynamics. Analysis includes:

Transmission rate β(t) and removal rate γ(t) are identified by solving using the inverse problem approach described in Section 4. They are calculated with P = 7, 10, and 14 day intervals over at least nine month period.

Reproduction rate R(t) = β(t)S(t)/γ(t), which represents the expected number of secondary infections produced by a single primary infected person. The effective reproduction rate R(t) is an important indicator of the fraction of the population that gets sick.

Simulated infectious cases I(t) using the A-SIR system. After estimating β(t) and γ(t) by solving the inverse problem, we find the numerical solution of the direct problem (1), (2), (3) with proper initial conditions. In all calculations, we use second order Runge–Kutta method, a type of predictor-corrector method. These simulations allow comparing the predicted values and the real data to test how well the method performs.

The time step τ is chosen to be τ = 0.01 for obtaining the numerical solution of the direct problem by Runge–Kutta method. The time-dependent parameters β, γ, and R are computed with τ = 1/40 and μ = τ. The experiments with different values of ϵ0 show that its optimal value is 10−8.

We chose the World and five Latin American countries in order to investigate the performance of the A-SIR model applied to an epidemic such as COVID-19. Brazil has the largest population and is the most affected country in Latin America; Argentina and Colombia are comparable in terms of population size although they have different climate and area; Dominican Republic and Honduras have smaller population size of approximately 10 million compared to Brazil, Argentina, and Colombia.

The simulations for COVID-19 epidemic are based on available data as of March 12, 2021, reported at (Worldometer, 2020). The reason we consider a period ending in March is because the vaccinations were relatively low in these countries until that date. It must be also noted that the available COVID-19 epidemic data does not include all infectious cases — this is outside the scope of this research.

The world

The dynamic of an epidemic in the World is important while studying this problem in a specific region. Fig. 4(a) presents the obtained values for the transmission rate β multiplied by the total population N, Fig. 4(b) – the removal rate γ, Fig. 4(c) – the reproduction rate R as functions of time, calculated over 7, 10, and 14 day intervals. The numerically obtained values for βN, γ, and R for the last day are given in Table 1. Fig. 4(d) shows the results for the infectives against the available data. The reproduction rate R had never been smaller than one until late January 2021. As seen in Fig. 4(d), the number of infected individuals worldwide was increasing until late January, when I(t) started to decrease.

Fig. 4

Estimated rates βN, γ, R, and the corresponding infectives for the World.

Table 1

Transmission rate βN, removal rate γ, and reproduction rate R for March 11, 2021 (World).

days	βN	γ	Re
7	0.018 825 6	0.018 514 7	1.001 296 4
10	0.019 172 6	0.018 851 5	1.001 537 7
14	0.018 489 5	0.018 762 8	0.970 419 4

Brazil

Fig. 5(a)–(d) show the numerically computed values of the time-dependent values of βN, γ, R, and I for Brazil. The values of βN, γ, and R for the last day, March 11, 2021, are given in Table 2. Brazil's reproduction rate has been increasing near the end of the period, being above one during the last few weeks.

Fig. 5

Estimated rates βN, γ, R, and the corresponding infectives (Brazil).

Table 2

Transmission rate βN, removal rate γ, and reproduction rate R for March 11, 2021 (Brazil).

days	βN	γ	Re
7	0.074 1	0.045 899 3	1.530 127 9
10	0.072 1	0.052 908 5	1.291 429 9
14	0.071 7	0.049 530 3	1.371 788 7

Argentina

Fig. 6(a)–(d) present the estimated values for the time-dependent βN, γ, R, and the corresponding infectives for Argentina, calculated over 7, 10, and 14 day intervals. Table 3 gives the specific values of βN, γ, and R on the last day of the period, March 11, 2021. Argentina's reproduction rate is close to one. This means that every sick individual is infecting one other individual.

Fig. 6

Estimated rates βN, γ, R, and the corresponding infectives (Argentina).

Table 3

Transmission rate βN, removal rate γ, and reproduction rate R for March 11, 2021 (Argentina).

days	βN	γ	Re
7	0.040 719 9	0.039 151 3	0.990 351 0
10	0.042 008 3	0.038 725 9	1.032 909 1
14	0.040 679 1	0.038 281 6	1.011 835 5

Colombia

The numerical results for the time-dependent parameters βN, γ, and R as functions of time for Columbia are presented in Fig. 7(a)–(c). The corresponding estimated infectives compared with the reported data are given in Fig. 7(d). Table 4 gives the values of βN, γ, and R for the last day of the considered period March 11, 2021. The effective reproduction rate oscillates around R = 0.85, meaning that the number of infectives I(t) is decreasing, as demonstrated in Fig. 7(d).

Fig. 7

Estimated rates βN, γ, R, and the corresponding infectives (Colombia).

Table 4

Transmission rate βN, removal rate γ, and reproduction rate R for March 11, 2021 (Colombia).

days	βN	γ	Re
7	0.087 596 6	0.100 194 7	0.835 023 0
10	0.088 240 7	0.100 988 6	0.834 550 5
14	0.087 084 7	0.094 022 3	0.884 641 3

Dominican Republic

The Dominican Republic reproduction rate R had been decreasing during the last weeks of the considered time period, as seen in Fig. 8(c). Fig. 8(a) and (b) display the obtained numerical values of the parameters β(t)N and γ(t), respectively. The transmission rate is decreasing while slightly oscillating. Table 5 gives the estimated values of βN, γ, and R for the last day of the considered period. Fig. 8(d) shows the number of infectives, also based on 7, 10, and 14 day intervals for approximating the rates β and γ.

Fig. 8

Estimated rates βN, γ, R, and the corresponding infectives (Dominican Republic).

Table 5

Transmission rate βN, removal rate γ, and reproduction rate R0 for March 11, 2021 (Dominican Republic).

days	βN	γ	Re
7	0.011 462 0	0.018 778 3	0.596 605 5
10	0.012 554 3	0.020 093 7	0.610 681 5
14	0.012 148 9	0.019 906 4	0.596 520 6

Honduras

The obtained numerical results for Honduras are presented in Fig. 9(a)–(d). The parameters βN and R are decreasing during the last weeks of the time period. The results based on 7, 10, and 14 day period are relatively consistent for Honduras. Table 6 presents the estimated values of βN, γ, and R for the last day of the 7, 10, 14 day time period. Clearly, although decreasing, the reproduction rate for Honduras is quite high in early March 2021, about 2.0. Fig. 9(d) shows the infectives I(t), again based on the last 7, 10, and 14 day time period used to estimate the transmission rate βN and the removal rate γ. Honduras has to further decrease their transmission rate in order to be able to slow down the COVID-19.

Fig. 9

Estimated rates βN, γ, R, and the corresponding infectives (Honduras).

Table 6

Transmission rate β, removal rate γ, and reproduction rate R for March 11, 2021 (Honduras).

days	βN	γ	Re
7	0.004 034 0	0.001 856 1	2.135 197 3
10	0.004 711 1	0.002 114 1	2.189 243 5
14	0.005 472 7	0.002 700 9	1.990 600 2

Comparison of real data and A-SIR and SIR predictions

This section presents simulations of the A-SIR model with time-varying coefficients and the classical SIR model with constant coefficients over a six week period (February 1, 2021–March 12, 2021), see Fig. 10, Fig. 11, Fig. 12, Fig. 13, Fig. 14, Fig. 15. The figures on the left display the results with the A-SIR model, while the figures on the right show projected values of infectives obtained by solving the classical SIR system. Both simulations solve the system (1)–(3) of equations using a second order Runge–Kutta method. The difference lies in the fact that the A-SIR model uses the estimated time-dependent rates, whereas the SIR model uses constant rates estimated as of February 1, 2021. The predicted values along with the real data for the same period, are presented in Fig. 10, Fig. 11, Fig. 12, Fig. 13, Fig. 14, Fig. 15 for the World, Brazil, Argentina, Colombia, Dominican Republic, and Honduras.

Fig. 10

Results from adaptive (left) and classical (right) SIR simulations compared with the real data for the last 41 days for the World.

Fig. 11

Results from time dependent (left) and classical (right) SIR simulations compared with the real data for the last 41 days for Brazil.

Fig. 12

Results from adaptive (left) and classical (right) SIR simulations and real data for the last 41 days for Argentina.

Fig. 13

Results from adaptive (left) and classical (right) SIR simulations compared with the real data for the last 41 days for Colombia.

Fig. 14

Results from adaptive (left) and classical (right) SIR simulations compared with the real data for the last 41 days for the Dominican Republic.

Fig. 15

Results from adaptive (left) and classical (right) SIR simulations compared with the real data for the last 41 days for Honduras.

The results show that the estimated coefficients in the A-SIR equations provide a reasonable description of the COVID-19 epidemic on the level of a country-wide distribution. Since the general trend is for a slow decrease of the effective reproduction rate R(t) in time, the real data is similar or slightly lower compared to the corresponding values obtained from the simulations.

The predictions for all countries are good. The reported values for Brazil oscillate while following the trend until February 22, 2021, and then they start to increase. There are oscillations in the reported data of Argentina, see Fig. 12. The A-SIR numerical results follow the trend of the data. The situation is similar with the results for Colombia, as shown in Fig. 13. The Dominican Republic experiences a slow-down in the reported infected cases during the last days of the six week period compared to the predicted values of the classical SIR model with coefficients estimated from February 1, Fig. 14. If the rates were taken from a later day (say, after February 7, 2021), then the predicted values would follow the real data. The A-SIR results for Honduras follow the trend closely, which the SIR projections are good at first and start to deviate after February 10, 2021 according to Fig. 15.

Conclusions

This work investigates the performance of a method for an epidemic based on an inverse problem approach for estimating the time-dependent transmission and removal rates in the A-SIR epidemic model. The inverse problem is solved by defining a minimization problem using the entire dataset for the examined population, with COVID-19 data available as of March 2021. The numerical results show that the rates used for finding the projected values of the infected cases using the adaptive SIR model are reasonably accurate. It is well-known that the classical SIR model possesses certain limitations in case of a long term infectious disease. Additionally, the assumption that the reported data include all the cases is a limitation, because it is not exactly true for COVID-19. Although the produced results are convincing, the match is not exact and the errors naturally increase over time. Small changes in the data (including errors) should not affect the results significantly.

Declaration of competing interest

No conflict of interests.