SIRV-D OPTIMAL CONTROL MODEL FOR COVID-19 PROPAGATION SCENARIOS Текст научной статьи по специальности «Математика»

EDN: OEDUSL УДК 51.77

SIRV-D Optimal Control Model for COVID-19 Propagation Scenarios

Viktoriya S. Petrakova* Vladimir V. ShayduroV

Institute of Computational Modelling SB RAS Krasnoyarsk, Russian Federation

Received 10.05.2022, received in revised form 26.05.2022, accepted 06.11.2022

Abstract. The article presents the compartmental differential formulation of SIR-type for modeling the dynamics of the incidence of viral infections, in particular COVID-19, taking into account the ongoing vaccination campaign and the possibility of losing immunity during some time period after vaccination or a disease. The proposed model is extended by considering the coefficients of the model as dependent on the social loyalty of the population to isolation and vaccination. This allows us to formulate the optimal control problem and build various scenarios for the development of the epidemiological situation. The results obtained on the basis of the considered models were compared with real statistical data on the incidence in the Krasnoyarsk Territory.

Keywords: scenarios of COVID-19 propagation, SIR-type model, optimal control model.

Citation: V.S. Petrakova, V.V. Shaydurov, SIRV-D Optimal Control Model for COVID-19 Propagation Scenarios, J. Sib. Fed. Univ. Math. Phys., 2023, 16(1), 87-97. EDN: OEDUSL

The COVID-19 pandemic has a huge impact not only on people's health, but also on the global economy and daily life. In response to this crisis, scientists, economists, and policymakers around the world are developing various countermeasures to contain the explosive spread of the virus, reducing the burden on the health system and the economy [1]. For example, measures such as travel bans, social distancing and masks have traditionally been chosen to contain the spread of infection. The invention of an effective vaccine significantly affected the dynamics of its spread.

Since the 1920s, mathematical models of the spread of infectious diseases have used as an important tool for analyzing epidemiological features and analyzing virus transmission [2]. There are a large number of different approaches, for example, more than 200 models have developed to predict the spread of COVID-19 for 2020-2021 [3]. Among them, the most commonly used models are compartmental ones of the SIR type (when the entire population is divided into susceptible to the virus, infected, and recovered people) and their modifications and adjustments [3]. For example, in [4], Caputo derivatives were used to formulate a fractional-order SIR-D (susceptible, infected, recovered, and dead) mathematical model for COVID-19 spread. Ramezani et al. in [5] proposed a modification of the SEIR-D (susceptible, exposed, infected, recovered, and dead) model to capture the behavior of the COVID-19 pandemic when asymptomatic infected individuals are taken into account.

However, as a comparison of predicted data with real ones shows, mathematical models of the spread of infectious diseases have a large number of sensitive parameters, the determination of which is a complex problem and is based on many assumptions. Unknown parameters must be estimated by fitting a model or solving inverse problems based on the data of the previous

* vika-svetlakova@yandex.ru https://orcid.org/0000-0003-1126-2148 t https://orcid.org/0000-0002-7883-5804 © Siberian Federal University. All rights reserved

time interval, that make forecast faulty [6]. Moreover, when modeling a disease for a long time, changing the parameters of the model is associated with many factors, for example, the introduction of quarantine restrictions, the mandatory wearing of masks, or an effective vaccination campaign.

This article proposes a modification of the standard SIR model by introducing into consideration the vaccinated part of the population and the possibility of the population to lose immunity after a while. In the second part of the article, the proposed model is extended to an optimal control problem, which made it possible not only to significantly improve the predictive properties of the model, but also to take into account social factors affecting the spread of the disease (compliance with coronavirus restrictions and vaccination) and build various scenarios for the development of the epidemic. The model was tested on statistical data for the Krasnoyarsk Territory for the spring of 2022.

1. Differential SIRV-D model

The most popular models for describing the spread of diseases are differential (compartmental) SIR models, which began in the work of the Scottish epidemiologists Kermack and McKendrick in the 1920s [2]. The essence of such models is the differentiation of the population into several parts (chambers) depending on the epidemiological status of individuals. Thus, in the simplest SIR [2] model, the entire population is divided into the following groups: "Susceptible" (S), the population that does not have immunity to the virus; Infected (I); and those who received immunity (R - "Recovered", in the simplest case, recovered and died individuals). The probabilities of an individual's transition from one chamber to another determine the system of differential equations described in many monographs and handbooks, for example, in [6]. Thus, the good analytical analyzability and computational simplicity of the SIR model led to the widespread use of this model, as well as its various modifications, which differ in the number of chambers into which the population is divided and the number of connections between them.

Some of the essential features of modeling the spread of coronavirus are the following statements:

• an individual does not receive permanent immunity and can be re-infected after weakening the immune system;

• at the time of model testing (2022), vaccination is the main antiviral measure;

• the dynamics of the disease varies depending on the region under consideration and is extremely unstable due to the emergence of new strains or the introduction of antiviral measures by the government;

• coronavirus has a number of symptomatic differences from other viruses, often (depending on the strain) is mild and can be asymptomatic.

The first and last two statements are specific features of the coronavirus infection and make it difficult to predict the dynamics of the epidemic for the medium-term (from 2 weeks to 2 months) and long-term (3 or more months) time periods using any known methods [3,6]. To take into account the first two features, we propose the SIRV-D model described by a system of differential equations

' dma

dt

dmi

dt

dmR

dt

dmv

dt

dm^

L dt

with initial condition

(t) = -Pma(t)mI(t) + $mR(t) — ema (t),

(t) = ¡3ma(t)mi(t) — jm,i(t) — Smi(t) + pmy (t) mi(t),

(t) = Y'mi(t) + umv (t) — 0mR(t),

■ (t) = —umv (t) + em s (t) — pmv (t) mi (t),

■ (t) = Smi (t)

mi(0) = m0i Wi G {S, I, R, V, D} .

(1)

(2)

The SIRV-D model (1), (2) is an extension of the simplest SIR model, where "dead" are separated from "immunized" and the group of "vaccinated" is introduced. A description of the model parameters is presented in Tab. 1, and a general flow diagram of the relationship between groups is shown in Fig. 1. Here, it is taken into account that the population loses its immunity and may become susceptible again, and the vaccinated part of population do not acquire immunity immediately, but after a certain time period. Note that yt G [0; T], where T is horizon time,

the law of conservation of the mass of the entire population ^ mi = 1 Vi G {S,I,R,V,D} is

i

satisfied.

Fig. 1. General flow of SIRV-D model

2. Estimation of model parameters

In order to correctly model the spread of the virus, in addition to an adequate basic model, it is necessary to correctly select its quantitative parameters. Some of these parameters are universal and the same in all territory, for example, the rate of loss of immunity or the rate of immunity formation after vaccination. The remaining parameters are specific to a particular region. Universal parameters are usually known and determined by the nature of the virus. Local parameters, as a rule, are obtained from the analysis of statistical data and the solution of inverse problems [8]. The task of determining the parameters is complicated due to the large amount of data, the variability epidemiological situation, the quality of collecting statistical data, determining the target functional (the question is what exactly should be interpolated: the number of infected people, the number of deaths, or the peak of incidence, etc.), as well the selection of solution algorithms emerging ill-posed problems. At the same time, due to the variability of the incidence dynamics, the determination of coefficients that approximate the real situation is possible only for a short time period (up to 2 weeks). More details about the difficulties encountered in determining epidemiological parameters and current methods for their

Table 1. The description of parameters of the SIRV-D model

Parameter description Symbol

The part of the population that is not immune at t time moment ms (t)

The part of the population that "active" infected at t time moment mI (t)

The part of the population that is immune at t time moment after suffering a mR(t)

disease or after being vaccinated

The vaccinated part of the population that has not yet received immunity at t mv (t)

time moment

Part of population who is died from the coronavirus by the time t mD (t)

Contagiousness ("infectiousness"), the probability of transmission of the virus in P

the population

The rate at which immunity is developed after infection Y

The rate of vaccination of the population £

Mortality a

The rate of loss of immunity, the probability of an individual becoming susceptible $

to the virus again

The rate of formation of immunity after vaccination w

The probability of vaccinated people getting sick before immunity develops P

solution can be found in [6,8-10]. Since in this work the SIRV-D model is only a start model for the optimal control problem presented below, and the modeling will be scenario-based, the accuracy in determining the coefficients of the model is not a priority task. So, to determine the coefficients of the SIRV-D model, we will assume that for each set of measurements on the t.-th day, the vector of unknowns qj = satisfies the following system of linear algebraic

equations Vi e {1,N — 1}

¡3lmlS mlI — YjmlI — 5lmlI = (m]+1 — mjfr)/2 — pmjV mlI,

YimiI = (ml+l — m/2 — wmV + $mR,

. . . . . . (3)

£j mSS = (mV+1 — m\- ^/2 + wmV + pmV m},

5imiI = (mi+1 — mi-1)/2 .

System (3) is obtained from the system of differential equations (1) after substituting the known measured values in t.-th day, where the central difference is taken instead of the first derivative. Note that instead of five equations from (1) to (3), only four equations are used, since otherwise, due to the properties of system (1), system (3) will consist of linearly dependent equations. Parameters w, $ and p will be considered known. Assume that the rate of loss of immunity $ is proportional to half a year, and the rate of formation of immunity after vaccination w is proportional to forty-two days. As the probability of infection for a vaccinated person without developed immunity, we take the data obtained by the Krasnoyarsk Medical University p= 18.3 % [11]. As parameters /3,j,£,S for further modeling, we choose the arithmetic mean over sets {Pj, ■ ■ -,Pj+k } , {Yj, ■ ■ ■, Yj+k } , {£j, ■ ■ ■, £j+k } , , ■■■, 5j+k }, where j > 1, j + k < N — 1. Tab. 2 shows the values of epidemiological parameters according to statistical data for the Krasnoyarsk city in the period from t1 =February 1, 2022 and tN-1 = 28 February 2022.

Fig. 2 shows a comparison of the approximation of statistical data for February in Krasnoyarsk with the solution of the SIRV-D model with constants from the first three rows of Tab. 2.

Note that any epidemiological model is interesting, as a rule, for predicting the development of epidemiological situation, that is, the model parameters are unknown for the entire predicted

Table 2. The value of the epidemiological constants for the Krasnoyarsk city (February, 2022)

Average for constants

ß Y e ö w P

01-28 Feb 0.22858 0.07385 0.00058 0.00045 0.00548 0.02381 0.18300

01-19 Feb 0.30832 0.08664 0.00057 0.00051

20-28 Feb 0.03222 0.04094 0.00064 0.00027

weekly average

first week 0.59188 0.20861 0.00009 0.00085 0.00548 0.02381 0.18300

second week 0.21143 0.02174 0.00068 0.00037

third week 0.08421 0.02278 0.00099 0.00029

fourth week 0.02679 0.04227 0.00059 0.00028

timeline, and can be selected as values from already known measurements for the previous period of time (sliding window method). The results of the predictive power of the SIRV-D model are shown in Fig. 3 (dots). So, to simulate the development of the epidemiological situation in the period from 02/07/2022 to 02/14/2022, the parameters obtained from the statistics of the previous week from 02/01/2022 to 02/06/2022, etc. were used. It can be seen that in this case, the SIRV-D model responds to a change in the situation with a delay and greatly overestimates the real situation. This property is inherent in any differential model of the SIR type.

3. Optimal control SIRV-D problem

As it can be seen from Tab. 2, the parameters 3, y are the most sensitive and determine the behavior of the curves that are the solution of model (1). Therefore, choosing them from the previous time interval for forecast leads to significant inaccuracies. In addition, the main criticism of SIR-type models is based on the assumption of population homogeneity without taking into account the differentiation of the population on social grounds: for example, age or compliance/noncompliance with antiviral restrictions. To overcome these difficulties, state the following control problem. Let aIs (t), aVc (t) are continuous functions denoting the cumulative loyalty of the population to isolation (observance of restrictions) and vaccination in a given region (isolation and vaccination strategies). Assume that \aI (t)|, \aV (t)| < 1 and the rate of spread of infection among the non-immune and vaccinated population 3, P and the rate of vaccination of the population e in (1) are now given as functions

3 = 3min + (3max " 3min) (1 + aIs (t)) /2;

P = Pmin + (Pmax - Pmin) (1 + O-Is (t)) /2; (4)

e = emin + (emax - emin) (1 + avc (t)) ¡2. We will find a set of functions ms, mI, mR, mv, mD, aIs, aVc that maximizes the functional

J = G(msiRVD,ais,avc)dt + $(msIRVD (0) ,mi (T)) (5)

■Jo

where the notation mSIRVD means the vector of state functions of the controlled system mSIRVD (t) = {ms (t), mI (t), mR (t), mV (t), mD (t)}. Here G determines the quality of the system functioning over the entire control interval, and $ are the terminal conditions, the final result of the control impact, determined by the combination of initial mi (0) and final mi (T).

Thus, we obtain the following optimal control problem: to maximize the functional (5) with constraints in the form of differential equations (1), taking into account (4). To solve the problem, we use the Pontryagin maximum principle. Choose smooth functions ^i (t) : [0, T] ^ R

Fig. 2. The comparison of statistical data with calculated one by the model (1). The gray line is statistical data, the dots present the SIRV-D solution obtained with average parameters for the entire simulated interval, the black line is the SIRV-D solution obtained with average parameters for two monotonous interval (from 1 to 19 February and from 20 to 28 February)

and multiply fa (t) Vi e {S, I, R,V, D} by the corresponding equation in (1). Write out the Pontryagin function [12]

H (msiRVD, aIs,aVc, ^sirvd) = AoG(msiRVD, aIs, avc)+ + ^s (t) (-ß (a i s) m s (t)mi (t) + fam,R(t) - £ (avc) ms (t)) + + ^i (t) (ß (ais) ms(t)mi(t) - ym-i(t) - 5mi(t) + p (ais) mv (t) mi (t)) + + ^R (t) (rmi(t) + umv (t) - fam,R(t)) +

+ 0V (t) (-umV (t) + e (aVc) ms (t) - p (aIs) mV (t) mI (t)) + 0d (t) (SmI(t)) .

Let {m*SIRVD (t), a}s (t), a'Vc (t)} be the optimal process in problem (1), (4), (5). Then there are simultaneously non-zero functions 0SIRVD (t) such that the following conditions are satisfied:

1. The Pontryagin function H (mSIRVD, aIs, aVc, 0SIRVD) of the problem for each t £ [0, T] reaches its maximum in aIs,aVc at the points a*Is (t), aVc (t) for mSIRVD = m*SIRVD (t),

0SIRVD = 0SIRVD (t) .

2. The vector function 0sirvd (t) satisfies the adjoint system of differential equations

30 = 3 (aI) mI(t) (0S (t) - 0I (t)) + e (aV) (0s (t) - 0v (t)) - dG/dmS , dt

dt = 3 (aI) ms(t) (0s (t) - 0i (t)) + y (0i (t) - 0r (t)) + +S (0i (t) - 0d (t)) + p (ai) mv(t) (0v (t) - 01 (t)) - dG/dmi , 30R ^ (0R (t) - 0s (t)) - dG/dmR ,

dt d0v dt d0D

p (ai) mi(t) (0v (t) - 0i (t)) + w (0v (t) - 0r (t)) - dG/dmv, -dG/dmD,

dt

with transversality conditions

^ d^jmj (0) ,mi (T)), ^ {T)= d^jmj (0) ,mt (T)) dmi (0) ' dmi (T)

To fulfill the first condition, we study (6) for the fulfillment of the necessary and sufficient conditions for the maximum of functions

dH dG + ßmax ., ßmin ms(t)mi(t) (0i (t) - 0s (t)) +

dais dais 2

+ Pmax - Pmin my (t)mi(t) (0i (t) - 0v (t)) = 0;

(9)

d2H d 2G , N

-2 =-2 < 0; (10)

d(ais) d(ais)

daH= da- + emx-mn w 0 w - *s w)=0 an -dLH-j = -^G, < 0. (12)

d(avc) d(avc)

From (10), (12) it can be seen that for the correct implementation of the maximum principle, the function G must be concave with respect to the control functions, then (9), (11) give the optimality conditions for the strategy chosen by the population.

4. Computational experiment

For the numerical implementation of the optimal control problem, we use the following iterative algorithm.

Input: initial values mSIRVD (0), initial control aIs (t) = aVc (t) = 0, set the maximum calculation precision eps.

Get a set of grid values mh(tk) by solving the system of differential equations (1) taking into account (4);

Get an approximate value of the cost functional Jh according to (5);

Until \J — max J| ^ eps do

(a) Get a set of grid values ^h(tk) by solving the system of differential equations (7) taking into account (4);

(b) Get a set of grid values ahs(tk),a^hc(tk) from equations (9), (11);

(c) Calculate the Pontryagin function (6) taking into account the current values mh(tk), fah(tk), ahs(tk), ahc(tk);

(d) Calculate the Pontryagin function (6) taking into account the current values mh(tk), '^i(tk), and boundary values ahIs(tk), a}hc(tk);

(e) As values ahs(tk), a!hc(tk) at the new iteration, select those that deliver a greater value to the function H;

(f) Get a new set of grid values mh(tk) by solving the system of differential equations (1) taking into account (4);

(j) Get an approximate value of the cost functional Jh according to (5);

Output: Select the grid functions mh(tk), ^h(tk), ahs(tk), aVc(tk) obtained at the last iteration as a solution of the optimization problem.

To solve systems of differential equations, the Runge-Kutta method of the 4-th order was used. To compare the predictive ability of the optimal control model and the SIRV-D differential setting on the data for the Krasnoyarsk Territory for February 2022, we consider a functional

J = j (d1(ms (t) + mI (t)+ mR(t) + mV (t))(1 + aIs (t)f —

Jo V (13)

—d2(1 + aVc(t))2ms(t) + ds(mI(t))2) dt,

where d1,d2,d3 are the weight constants and £,p are the constants vary within £ e [O^OOOBQ, 0^00099], p e [0.1464, 0.2196].

Let's consider several scenarios.

1. Optimistic scenario. We minimize the functional (13), where the epidemiological parameter P e \pfifit, fit] and Pt is table value used in the forecast interval.

2. Pessimistic scenario. We maximize the functional (13), where P e [Pt, 1Apt].

3. Average scenario. We minimize the functional of the form

J = J ('d,1(ms (t) + mI (t) + mR(t) + mV (t))(aIs(t))2 + d2(aVc(t))2ms (t)+ (14) +P(aIs(t))ms(t)mI(t) + SmD(t) — £(aVc(t))mV(t)) dt, where parameters P e [06pt, 1Apt].

Minimization of the functional (13) leads to the fact that it is beneficial for the population to be loyal to isolation and vaccination, which leads to a decrease in the number of infected people. The pessimistic scenario leads to the opposite. The average scenario (14) leads to the selection of isolation and vaccination strategies by the population to minimize the daily increase in cases and deaths.

Fig. 3 shows a comparison of the predictions of a simple differential SIRV-D model and an optimal control model based on it. As for the SIRV-D model, the model parameters for the

Fig. 3. The comparison of statistical data with forecasting of simple differential model SIRV-D (1) and optimal control SIRV-D (1), (4), (5), (7), (8). For the first week of March the prognose for different scenarios is presented

forecast for the new week were chosen as the values obtained from the previous week. The scenario branching is shown for the first week of March 2022. The black curve is made up of the scenarios that showed the best approximation with the statistical data. It can be seen from the figure that the optimal control model can give a more accurate forecast, but it requires an adequate assessment of the situation when constructing a new forecast scenario.

However, one of the drawbacks of the proposed optimal control model is the fact that not only the sensitive parameters 3, P, which are responsible for the probability of becoming infected, are corrected, but also the parameter e, which varies within small limits. While the population recovery rate parameter y, which has a significant impact on the dynamics of the spread of infection cannot be corrected by considering it as dependent on the control function. One way to solve this problem would be to consider the basic reproductive number of the population [2,3]

Ro = P/j .

(15)

It is obvious that for simple dependences if R0 > 1, then there is an increase in the incidence, otherwise there is a decrease. In order to build more accurate forecasts, one can consider a change ¡3 within [7 — ¿1; 7 — S2]. Fig. 4 shows a forecast of a new increase in the incidence in the Krasnoyarsk Territory in 2022 with a pessimistic forecast (since at May of 2022, all coronavirus measures and mandatory vaccination have been canceled), provided that 3 € [0.4y; 1.67].

Fig. 4. New wave prognoses by model (1), (4), (5), (7), (8)

The model shows that the next significant increase in incidence may occur by middle of July this year. However, as already shown above, a correct long-term forecast is possible only if the epidemiological situation is preserved (no new strains or restrictions, as well as maintaining the mobility of the region's population).

5. Conclusions

As a result, two well-analyzed analytical models were obtained to describe the possible dynamics of the disease, taking into account vaccination campaign. On the basis of a simple differential model, an optimal control model is constructed, where some coefficients are dependent on the social moods of the population. It is shown that, unlike the differential model, the optimization problem can give a more accurate forecast, and also allows you to quickly take into account the introduction of restrictive measures by changing the cost functional. Several scenarios of the dynamics of the epidemiological situation for the Krasnoyarsk Territory have been constructed, depending on the different social behavior of its inhabitants.

The work was supported by the Russian Science Foundation grant no. 20-61-46017.

References

[1] D.Rafiq, A.Batool, M Bazaz, Three months of COVID-19: a systematic review and metaanalysis, Rev. Med. Virol., 30(2020), no. 4, e2113. DOI: 10.1002/rmv.2113

[2] W.O.Kermack, A.G.McKendrick, A contribution to the mathematical theory of epidemics, Proc. R. Soc. Lond. Ser. A, Contain. Pap. Math. Phys. Character, 115(1927), 700-721.

[3] O.I.Krivorotko, S.I.Kabanikhin, Mathematical models of COVID-19 spread, Preprint: https://arxiv.org/pdf/2112.05315 (December 2021), 700-721.

[4] H.Mohammadi, S.Rezapour, A.Jajarmi, On the Fractional SIRD Mathematical Model and Control for the Transmission of COVID-19: the First and the Second Waves of the Disease in Iran and Japan, ISA Transactions, 124(2021), 103-114. DOI: 10.1016/j.isatra.2021.04.012

[5] S.B.Ramezani, A Amirlatifi, S.Rahimi, A novel compartmental model to capture the nonlinear trend of COVID-19, Comput. Biol. Med., 134(2021), 104421. DOI:10.1016/j.compbiomed.2021.1044

[6] O.I.Krivorotko, S.I.Kabanikhin, N.Y.Zyatkov et al., Mathematical modeling and forecasting of COVID-19 in Moscow and Novosibirsk region, Numer. Analys. Appl, 13(2020), 332-348.

[7] F.Brauer, P.Driessche, J.Wu, Mathematical Epidemiology, Lecture Notes in Mathematics, 355(2008), 412.

[8] T.T.Marinov, R.S.Marinova, Dynamics of COVID-19 using inverse problem for coefficient identification in SIR epidemic models, Chaos, Solitons and Fractals, 5(2020), 100041. DOI:10.1016/j.csfx.2020.100041

[9] A.Cintron-Arias, H.T.Banks, A.Capaldi, A.L.Lloyd, A sensitivity matrix based methodology for inverse problem formulation, Journal of Inverse and Ill-Posed Problems, 17 (2009), no. 6, 1-20. DOI: 10.1515/JIIP.2009.034

[10] A.Paticchio, T. Scarlatti, et al., Semi-supervised Neural Networks solve an inverse problem for modeling Covid-19 spread, Preprint: arXiv:2010.05074v1 [cs.LG], 2020, 6.

[11] Online: https://ngs24.ru/text/health/2021/09/23/70151987/ (Last access - 05/21/2022).

[12] A.D.Lewis, The Maximum Principle of Pontryaginin control and in optimal control, Interim Reportson work of the International Institute for Applied Systems Analysis, Vol. 143, 2006.

Модель оптимального управления SIRV-D для построения сценариев распространения COVID-19

Виктория С. Петракова Владимир В. Шайдуров

Институт вычислительного моделирования СО РАН Красноярск, Российская Федерация

Аннотация. В статье представлена камерная дифференциальная постановка SIR-типа для моделирования динамики заболеваемости вирусными инфекциями, в частности COVID-19, с учетом проводимой кампании вакцинации и возможности потери иммунитета через некоторый промежуток времени после вакцинации или заболевания. Предлагаемая модель расширена за счет учета коэффициентов модели как зависимых от социальной лояльности населения к изоляции и вакцинации. Это позволяет сформулировать задачу оптимального управления и построить различные сценарии развития эпидемиологической ситуации. Результаты, полученные на основе рассмотренных моделей, сравнивались с реальными статистическими данными о заболеваемости в Красноярском крае.

Ключевые слова: сценарии распространения COVID-19, модель SIR, модель оптимального управления.