OPTIMAL CONTROL IN THE NETWORK MODEL OF BI-VIRUS PROPAGATION Текст научной статьи по специальности «Математика»

Contributions to Game Theory and Management, XV, 265—286

Optimal Control in the Network Model of Bi-virus

Propagation *

Liu Xiuxiu and Elena Gubar

St. Petersburg State University, Faculty of Applied Mathematics and Control Processes, 7/9, Universitetskaya nab., St. Petersburg, 199034, Russia E-mail: st091520@student.spbu.ru; e.gubar@spbu.ru

Abstract This thesis is devoted to the spread of virus in human society. First, we modified the original basic epidemiological model, and divide the system into M different types. The optimal control problem is formulated for the changes of compartments in different states. The total cost is then minimized, the Pontryagin maximum principle is used to solve this nonlinear optimal control problem. Next, we prove that the optimal policy has the simple structure. Finally, we fit the propagation process of this model using Matlab. Consider the situation where there are only two types of virus in the system, and compare the two types of virus appear at the same time and at different times.

Keywords: optimal control, virus propagation, epidemic model. 1. Introduction

The basic models of epidemics are the susceptible-infected model (SIS) and the susceptible-infected-recovered model (SIR) (Bichara et al., 2014). The total population is divided into a susceptible, an infected, and a recovered subpopulation. We assume that the general progression of an epidemic is a dynamic process, moving from the susceptible to the infected and finally to the recovered population.

Viruses are neither living nor abiotic, but parasitic self-replicator. They can infect all living organisms with cellular structures. When a susceptible compartment encounters an infected compartment, there is the potential for virus to spread.

On the one hand, compartments throughout the system are divided into many different types, virus of different types can spread randomly throughout the system. On the other hand, we assume that virus of different types appear at different times, and that vulnerable compartments are likely to be infected by different virus at different times (Gubar and Zhu, 2013).

Mathematical models based on nonlinear differential equations have been developed and applied to various systems, and we can model all systems with epi-demiological behavior (Khouzani et al., 2011). First, an optimal control formula is established according to the changes of compartments in different states. The total cost is then minimized, resulting in an ideal tradeoff between network security and resource consumption. Finally, the Pontryagin maximum principle is used to solve this nonlinear optimal control problem.

*This work was supported by the China National Scholarship for Study Abroad. https://doi.org/10.21638/11701/spbu31.2022.20

2. System Model 2.1. Dynamics

We first determine the state evolution and then formulate the system model. A system consists of N compartments, and at time t, the number of susceptible, exposed, infected, recovered, and dead states is ns(t),nE(t),nI(t),nR(t) and nD (t), the corresponding fractions are S(t) = ns (t)/N,E(t) = nE (t)/N,I(t) = ni (t)/N,R(t) = nr(t)/N,D(t) = nD (t)/N. So for all t, we have S(t) + E(t) +1(t) + R(t) + D(t) = 1. For type i, at time t, the fractions for the states susceptible, exposed, infected, recovered, and dead are Si(t),Ei(t),Ii(t),Ri(t),Di(t), respectively. In our model, we assume that the compartments of each type are stable and do not change over time.

Susceptible compartments are compartments that can be easily attacked by viruses but have not been infected; exposed compartments are compartments that have been infected but do not have the ability to spread; infected compartments are compartments that have been contaminated by virus; recovery compartments are compartments that have become immune to virus; dead compartments are lifeless (Eshghi et al., 2016).

The compartments of the system can be stratified into different types of M (Eshghi et al., 2016). The number of compartments in these types need not be exactly equal. Compartments of type i contact compartments of type j at the rate Hji.

In each type, a set of compartments filled with the vaccine, called the immune compartment group, is preselected. The immune swarm can distribute vaccines to vulnerable compartments and infected compartments, immunizing vulnerable compartments against the virus and potentially curing infected compartments. When a susceptible compartment becomes immune or an infected compartment recovers, it is converted to a recovery compartment. The number of immune individuals is NR, where 0 < R < 1. We assume that immune individuals are not infected, so they regain their state from the start. At t = 0, let 0 < S(0) = So < 1,0 < E(0) = Eo < 1,0 < I(0) = Io < 1,D(0) = Do = 0. Thus, R(0) = 1 — S0 — E0 — I0 — D0. If no virus exists, I0 = 0.

Fig. 1. State transitions

In modeling, we assume that susceptible compartments that have touched exposed compartments are immediately exposed and homogeneously mixed, i.e., exposed compartments have the same probability of being in contact with any other compartment. The susceptible compartment of type i to contact the exposed com-

partment of type j at the rate j. Let an exposed compartment of type i contact an infected compartment of type i at a rate w,. Immune compartment group distribute their vaccine on contact to infected and susceptible compartments, with i type immune compartment group distributing vaccine at u,. Obviously, 0 < u, < 1. Immune compartment group can touch susceptible or infected compartments at a rate of W. The mortality rate of each infected compartment of type i is S,, here, the mortality in the other compartments are negligible. Where 0 < j < 1,0 < w, < 1,0 < W < 1 and 0 < Si < 1.

The efficacy of the vaccine on infected compartments may be less than on susceptible compartments. If the vaccine heals the infected compartment, the infected compartment's state transitions to recovered, otherwise it remains infected. 9, represents the effectiveness of patching i type of infected compartment.

Figure 1 illustrates the transitions between different states of compartments and the notations used in type i.

If the total number of compartments (N) is large, then S(t), E(t), I(t), R(t) and D(t) converge to the solution of the following system of differential equations:

M

Si = — y ] Sj^jjEj — SiWjRjUi j=1

M

Ej = Sj^jjEj — EjWjlj (1)

Ij =EjWjIi — IjWiRiUi 9, — Sjlj

Rj =SjWjRiUi + IjWiRiUj9i,

where, by writing S(t) = YM= i Sj(t), E(t) = EM=i Ej(t), I(t) = £M=i I(), D(t) = M=1 Di(t), R(t) = Y, = Ri(t), the initial conditions and state constraints

are given by

0 < S(0) = S0 < 1,0 < E(0) = E0 < 1,0 < I(0) = I0 < 1, D(0) = Do = 0, R(0) = 1 - So - Eo - Io - Do, S(t) > 0, E(t) > 0, I(t) > 0, R(t) > 0, D(t) > 0, S(t) + E(t) +1 (t) + R(t) + D(t) = 1.

(2)

2.2. The Objective Function

We seek the minimum total cost of the system. Since the malicious activity of the virus affects the evolution of the system, the system incurs a cost of f (I(t)), g(E(t)), k(D(t)) at each time t. The network also benefits at a rate of L(R(t)) due to the removal of uncertainty about the state of the vaccine. where f (-),g(-),k(-),L(-) are all non-decreasing and differentiable functions such that f(0) = g(0) = k(0) = L(0) = 0. Note that it is natural to assume that f (.),g(.), k(.) can be any function.

Suppose that each active immune compartment group at time t consumes resources at a rate of h(u(t)), because it uses a rate of u(t) at time t. Here h(0) = 0 and h(x) > 0 if x > 0. Suppose the cost of consuming additional resources at time t is given by a sum of the form ^M=1 Rihi(ui).

Di =àiii

Using the above parameters, the total cost is given by an expression of the form

T m

J = (f (I) + g(E) + k(D) - L(R) + Y Rihi(ui))dt. (3)

0 i=i

3. Optimal Protecting Strategies

3.1. Theoretical for computing the optimal controls

We apply Pontryagin's Maximum Principle to obtain a framework for solving the optimal control problem as posed (Grass et al., 2008). Let ((S, E, I, D, R),u) be an optimal solution to the problem posed in the previous section, consider the Hamiltonian H, and the corresponding co-state or adjoint functions AS, Af, A1, Af and AR, defined as follows:

M

H(u(t)) =f (I(t)) + g(E(t)) + k(D(t)) - L(R(t)) + Y Rihi(ui)

i=i

M

+ Y(AS Si + Af Ei + Af Ii + Af Di + aRR i)

i=i

M

=f (I(t)) + g(E(t)) + k(D(t)) - L(R(t)) + Y Rihi(ui)

i=i

M M

+ Y((-AS + Af) Y SiHiEj + (-AS + AR)SiOJiRiUi

i=1 j=1

+ (-Af + Af )EiUiIi + (-Af + AR )IiUJiRiUidi + (-Af + Af )Si Ii),

Where the adjoint functions AS, Af, A1, Af and AR are continuous functions, We

have a differential equation,

dH M

ASS = -^ = (ASS - Af )Y VnEi + (ASS - ARHRiU

i j=i

Af = - H = - dj§E1 + Y(AS - Af )mSj + (Af - Af ).iIi i i j=i dH df (T) (5)

Af = - Wi = - f1 + (Af - A1 )EiUi + (A1 - AR)cJiRiUiOi + (Af - Af )Si

■ D dH dk(D) A D

(4)

dDi BDi dH dL(R)

AR = --R = -LR - hi(ui) + (AS1 - AR)SiCJiUi + (Al - AR)IicJiUiûi,

Together with the final conditions

Af (T) = 0, Af (T) = 0, A'(T) = 0, Af (T) = 0, AR(T) = 0. (6)

Then PMP implies that the optimal control at time t satisfies the following conditions.

ui G argminH(v), (7)

where the minimization is over the space of admissible controls.

From (3), it can be shown that vector minimization can be expressed as scalar minimization

Ui(t) € argminRi(t)Yi(x,t), (8)

0< x< 1

where,

Yi(x,t) = hi(x) - &(t)x, (9)

&(t) = (Af - AR)SnJi + (Af - *?)IiuJiOi. (10)

First, from (6) that Yi(0,t) = 0, hence (8) results in Ri(t)ii(ui,t) < 0. That is, for all t,

Yi(ui,t) < 0. (11)

3.2. Structure

We are now ready to identify the structures of the optimal controls (ui(t),..., um (t)):

Lemma 1. For any i, ^(t) is a decreasing function of t.1

In economic terms, the adjoint function represents the shadow price. In simple terms, Af should be a positive number because it represents the additional cost per unit time that the system incurs as the proportion of dead compartments increases. Similarly, Af should also be a positive number, and Af should also be a negative number. Therefore, we expect Af - Af > 0,Af - Af > 0.

Lemma 2. The positivity constraints Af (t) - Af (t) > 0 and Af (t) - Af (t) > 0 hold for all i =1,...,M and all t e [0, T) .2

Theorem. Assuming the existence of an optimal control, for types i such that R > 0: if hi(■) is concave, then the optimal control for type i has the following structure: u*(t) = 1 for 0 < t <ti, and u*(t) = 0 for ti < t < T, where ti e (0, T). If hi(-) is strictly convex, then the optimal control for type i, ui(t) is continuous and has the following structure: u*(t) = 1 for 0 < t < tfu *(t) = 0 for t\ < t < T, and u*(t) strictly decreases in the interval (t1,t2), where 0 < t1 < t2 < T (Gubar and Zhu, 2013, Eshghi et al, 2016).

Proof. This time the Hamiltonian

M M

\S , \E\

H(u) = f (I) + g(E) + k(D) - L(R) + £((-A? + XE) E SiHiEj

i=1 j=1 M

+(-XE + XI)EitViIi + (-XI + Xf )SiIi) + Y^(Ri(hi(ui) - tiUi)).

i=1

(12)

1) Let hi(■) be a concave function, i.e., hi (■) < 0, then the Hamiltonian is a concave function of ui,i = 1,2,...,M. There are two different possibilities for u.i e [0,1] that minimize the Hamiltonian(Gubar et al., 2021),

1See Appendix 1 for proof of Lemma 1.

2 See Appendix 2 for proof of Lemma 2.

If at the time t

H(0) > H(1)

that is, 0i > hi,(1), then the optimal control is u* = 1; otherwise u* = 0.

0i (t) is a decreasing function of t, and Lemma 1 will prove this. Define 0. (ti) = hi (1), We have

' 0i(t) >hi(1),t e [0, ti)

0i(t) <hi(1),t e [ti, T]. (13)

So,

'u* = 1,t e [0,ti) (14)

u* = 0,t e [ti,T]. (14)

2) Let hi(•) be a strictly convex function, i.e., hi (•) > 0, then the Hamiltonian is a convex function of ui ,i = 1, 2,...,M. There are three different possibilities for ui e [0,1] that minimize the Hamiltonian,

If at time t,

Fig. 3. The case of convex function - h (•) > 0 d(hi (u.) - 0iU.)

dui

then the optimal control is u* = 0.

If i

d(hi (u.) - 0.u.)

dui

then the optimal control is u* = 1. Otherwise,

d(hi (u.) - 0.u.)

k=0 = hi (0) - 0i > 0,

|ui=i = hi(1) - 0i < 0,

\ui =u* — 0

dui

we can find such a value u* e (0,1).

0. (t) is a decreasing function of t. Define 0. (ti) = hi (1), 0. (t?) = hi (0). We have

' 0.(t) > h'%(1),t e [0, ti]

h.(0) <0.(t) <h.(1),t e (ti,t?) (15)

0i(t) < hi(0),t e [t?,T].

So,

u* = i,t e [0,t}]

u* = h— fa),t e (t},t\)

u* = 0,t e №,T].

(16)

4. System Model of Two Types

In this section, we consider a case of M = 2 and it is an extension of the model considered in the previous section. In order to better and more intuitively describe the changes of the system state when different viruses appear at different times, only the case of M = 2 is considered here. The spread of viruses is a process from local to global. The case of M = 2 means that at a certain time, a particular virus is found in two different local environments, and the spread of these two different viruses will affect the overall process. Compartments of type i contact compartments of type j at the rate of j. In different types, the totals may not be equal, and of course the value of each parameter may vary.

1

2

Fig. 4. State transitions of two types

Figure 4 illustrates the thansition schemes for SEIRD model with two types of virus. If the total number of compartments (N) is large, then fractions S (t), E (t),I (t), R(t)

and D(t) converge to the solution of the following system of differential equations:

51 = Ei — S 1^21 E2 — S1UJ1R1U1

52 = —S2^12 Ei — S2^22E2 — S2UJÍ R2U2

Ei = Si^iiEi + S2 ^12 Ei — Ei^iIi

E2 = Si^2iE2 + S2 ^22E2 — E2U2I2

11 = Eiwi Ii — IiüiiRiuiOi — Sili

12 = E2W212 — I2W2R2U2O2 — S2I2 Di = SiIi

D2 = S2I2

Hi = SiWiRiui + IiujiRiuidi R2 = S2W2R2U2 + I2W2R2U2O2

(17)

where, by writing S(t) = S1(t) + S2(t),E(t) = E1(t) + E2(t),I(t) = I1(t) + h(t),D(t) = Di(t) + D2(t),R(t) = Ri(t) + R2(t), the initial conditions and state constraints are given by

(18)

0 < S(0) = S0 < 1,0 < E(0) = E0 < 1,0 < I(0) = I0 < 1, D(0) = Do = 0, R(0) = 1 - So - Eo - Io - Do S(t) > 0, E(t) > 0, I(t) > 0, R(t) > 0, D(t) > 0, S (t) + E (t) + I (t) + R(t) + D(t) = 1.

The aggregate cost is given by an expression of the form:

T

J = J (f (I) + g(E ) + k(D) — L(R) + Rihi(ui) + RJ2h2(u2))dt. (19)

We apply Pontryagin's Maximum Principle to obtain a framework for solving the optimal control problem as posed. Consider the Hamiltonian H, and corresponding co-state or adjoint functions AS, Af, A1, Af and AR, defined as follows:

H(u(t)) =f (I(t)) + g(E(t)) + k(D(t)) - L(R(t)) + Rihi(ui) + R2h2(u2)

+ (-Af + Af )(Si mEi + S\mE2) + (-Af + A?)Si<JiRiui

+ (-Af + A[)Ei uili + (-A{ + AR )InJiRiuidi + (-A{ + Af )Si Ii

+ (-Af + Af )(S2 Mi2Ei + S2V22E2) + (-Af + AR)S20J2R2u2

+ (-Af + Ai)E2 U2I2 + (-A2 + AR)I2 Ü2R2u202 + (-A2 + Af )Ö2h,

(20)

We got differential equation,

A? = (A? - Af )(MIIE1 + M2iE2) + (A? - X1R)UJ1R1U1 A? = (A? - Af)(M12EI + ^22^2) + (A? - AR)CJ2R2U2

Af = -dE) + (A? - AE)MiiSi + (A? - Af)Mi2^2 + (Af - A{)uJi AE = -^^ + (A? - Af )M2iSi + (A? - Af )M22^2 + (Af - A2)^2/2

0E2

= -^ + (Af - Al)Ei^i + (Al - AR)cJiRiui6i + (A{ - Af

Oil

\ I

A ?

\ ? — A2 — -

^ffi + (Af - a2 )E2^2 + (a2 - \R)uJ2R2U2e2 + (a2 - a?? )Ö2 (21)

012 dk(D) dD1 dk(D)

dD■

2

AR = ^LR^ - hi(ui) + (A? - All)Siumi + (Al - AR)hu>iuiQi = ^LRR - h2(u2) + (A? - AE)S2UJ2U2 + (A2 - AR)l2U>2U292,

along with the final conditions

A?(T) = 0, Af (T) = 0, A'(T) = 0, Af (T) = 0, AR(T) = 0 (22)

Then PMP implies that the optimal control at time t satisfies

ui € argminH(vi) U2 € argminH(v2).

(23)

where the minimization is over the space of admissible controls.

We modified the base model, and on the modified model, we considered that the system contains multiple types with different subpopulations. Our main work is to study the propagation behavior of virus when the system contains only two subpopulations, here, this is an uncontrolled model.

5. Numerical Investigations

In this section, we confirm our results with some numerical simulations. For the theoretical model introduced in the previous section, it is an uncontrulled model. For simplicity and without loss of generality, we only simulate the case where the cost function h(u{) is concave. The focus of this section is to study the behavior of the modified model, we will conduct two experiments: the first, comparing the presence of two different virus at the same time and at different times, on the dynamics of susceptible, exposed, infected, recovered, and dead compartments to evaluate. The second experiment simulates the effect of different parameter values on virus propagation using the first set of data from the first experiment.

2

Fig. 5. State transition of only one type of in the system

First, we assume that there is only one type of subpopulation in the system, only one type is infected, and the state transition is shown in Figure 5.

Let Ii = 0.0001, Ei = 0.00005, Si = 0.99, Di = 0, for i = 1, also, = 0.76, wi = 0.55, cJi = 0.43; ui = 0.55, 0i = 0.65, ¿i = 0.01.

Fig. 6. The changing trend of each state for only one type of virus in the system

According to Figure 6, the susceptible fraction decreases rapidly when t e [0, 28], and then gradually flattens. The exposed fraction grow rapidly at t e [0, 28], and reach the highest point at t = 28. At this time, the fraction of exposed compartments is 0.6746. And then decreases rapidly, and gradually flattens after t = 67. The infected fraction has no obvious change when t e [0, 28], and it grows rapidly when t e [28, 56], and reaches the highest point when t = 56. At this time, the fraction of the infected compartment is 0.3481. After that, it decreased rapidly, and gradually flattened after t = 86. The dead fraction does not change significantly when t e [0, 45], and then slowly rises. After t = 117, it is basically in a stable state. The recovery fraction increases significantly when t e [0, 25] U [42, 95], while the increase in t e [25, 42] is small, and after t = 95, it is basically stable.

In the early stage of virus propagation, most compartments are in a susceptible state, so the decline rate of susceptible fraction is relatively fast, and these compartments just turn into exposed compartments or recovered compartments. Therefore, exposed and recovered fraction grow rapidly at this stage. When the exposed compartments grow to a certain number, some of them are converted into infected compartments, and the infected fraction increase significantly. The growth of exposed and infected fraction is bound to reach a commanding height at some point, when virus spreads the most. After the commanding heights, the proportion of exposed and infected fraction will drop significantly. As the end compartments of the system, the death and the recovery fraction will slowly increase in propor-

SEIDR model

50 100 1 50 200 250 300 35Ö 400 Time

Llir.QJMHS)

SEIDR model- type 1 SEIDR model - type 2

0 50 100 150 200 250 300 350 400

c)

Fig. 7. Same time case of the first data: a) only type 1 in the system; b) only type 2 in the system; c) type 1 and type 2 in the system

tion as the virus spreads, and reach a stable state after a certain time. Recovery fraction grow slowly when t E [25,42], because almost all susceptible compartments have been converted to exposed compartments or recovered compartments at this time, while infected fraction have not grown significantly. The proportion of susceptible and infected fraction is small, and only a small number of compartments can be converted into recovery compartments. Over time, compartments transition between different states, and the number changes significantly. In the late stage of virus propagation, the number of compartments in all states tends to stabilize.

Then we discuss the case where the system contains two different types. In the two distinct types of subpopulations, state transitions not only occurred within, but were also affected by mutual infection between them. In different types, the totals may not be equal, and of course the value of each parameter may vary. At this point, we consider two different scenarios, where two types of virus appear at the same time, and the two types of virus appear at different times. The state transition is shown in Figure 4.

For the first case, we assume that the virus all appears when t = 0, and we can fit the corresponding change graph. Let Ni = 50000, N2 = 100000, N = Ni + N2,Ii = 0.0002, I2 = 0.00015, Si = 0.998, S2 = 0.999, Di = 0; D2 = 0,Ei = 0.00002, E2 = 0.00002, ^ii = 0.76, ^i = 0.35, ^22 = 0.56, ^ = 0.35, ¿i = 0.55,(2 = 0.55, (¿1 = 0.43,^2 = 0.45, ui = 0.55, U2 = 0.45,0i = 0.65,02 = 0.4, Si = 0.01, ¿2 = 0.015. Shown it in Figure 7.

Then we try to change the value of each parameter and study the same characteristics of the transitions between states when the parameter values are different.

SEIDR model - type 1 SEIDR model - type 2

0 100 200 300 400 500 600 700 8C 0 100 200 300 400 500 600 700 800

Time Time

a) b)

SEIDR model-type 1 8,2

1 f> -nft n QOjl^-1-1-1-1-

0 100 200 300 400 500 600 700 800 Time

c)

Fig. 8. Same time case of the second data: a) only type 1; b) only type 2; c) type 1 and type 2 in the system

Let N1 = 200000, N2 = 100000, N = N1 + N2,/1 = 0.00005,12 = 0.0001, Si = 0.9995, S2 = 0.999, Di =0; D2 = 0, Ei = 0.000025, £2 = 0.00002, = 0.56, ^21 = 0.45, ^22 = 0.66,^12 = 0.55, = 0.55,^2 = 0.74, CJ1 = 0.42, CJ2 = 0.36, « = 0.4, «2 = 0.3,01 = 0.7,6>2 = 0.4, ¿1 = 0.03, ¿2 = 0.01. Shown it in Figure 8.

According to a) and b) of Figure 7 or Figure 8, we find that when only type 1 or type 2 is included in the system, the transitions between different states are roughly similar to those in Figure 6. Since the parameter values selected by the two are not equal, the steady state values of the corresponding peaks are not equal. We mainly discuss c) of Figure 7 or Figure 8, which shows that when the system includes two different types, virus not only spreads within the same type, but also interacts between different types, and virus from other types will susceptible compartments are converted to other types of exposed compartments.

Table 1. Fraction of each state at steady state

R D

Only type 1 0.8514 0.1484

Only type 2 0.5922 0.4078

Type 1 and Type 2 0.3158 0.6831

According to Figure 9 and Table 1, when there is only one type in the system, the recovery fraction is higher than the death fraction at steady state. However,

SEIDR model SEIDR model

a) b)

Fig. 9. Same time case of the first data: a) recovered fraction; b) death fraction

when there are two types in the system, the junctions between the two types also interact, resulting in higher death fraction than recovery fraction.

Table 2. Fraction of each state at the peak

data 1 data 2

E I E I

Only type 1 Only type 2 Type 1 and Type 2 0.9025 0.7218 0.9251 0.7892 0.9738 0.5725/0.5679 0.981 0.7852 0.9883 0.9354 0.9943 0.6318

According to Figure 10 and Table 2, after the interaction of two types of compartments in the system, the fraction of exposed compartments at the peak are slightly higher than those of the system with only one type, and the fraction of infected compartments at the peak are slightly lower than those of the system with only one type type of situation. Because when there are two types in the system, the exposed compartments come from two types of susceptible compartments, not just from the same type of susceptible compartments under the action of different virus from the two types. The ^-values of the two types are different, causing the infection compartments to peak at t = 67 and t = 92, respectively.

Table 3. Time required for each state to stabilize

data 1 data 2

S I E D R S I E D R

Only type 1 Only type 2 Type 1 and Type 2 28 100 50 102 114 37 150 50 196 238 42 300 100 343 361 34 159 63 194 203 30 369 45 453 485 35 489 159 587 536

According to Table 3, when there are two types in the system, the time required for each state to stabilize is longer than when there is only one type in the system. The problem we face is exacerbated when two types of virus are present in the system.

Fig. 10. Same time case of the first data: a) exposed fraction; b) infected fraction. Same time case of the second data; c) exposed fraction; d) infected fraction

Then, we fit the occurrences of the two types of virus at different times. This situation is more realistic, and it is very unlikely that two different pieces of virus will appear at the same time. We assume that for type 1, the virus appears at t = 0, and for type 2, it appears at t = 180.

According to a) and b) of Figure 11, virus of type 1 appears at t = 0, and each state changes significantly when t G [0,100]. Virus of type 2 appears at t = 180, and each state changes significantly when t G [180,350]. c) of Figure 11 Shows the changes in the fraction of compartments incompartment each state when two types of virus appear in the system at different times. When type 1 virus appears, it will cause the susceptible fraction to decrease at t G [0, 76], the exposed fraction to increase at t G [0, 70], to decrease at t G [70,144], and the infected fraction to increase at t G [70,125], decreasing at t G [125, 208]. The susceptible fraction remain stable in t G [76, 176], because the number of susceptible compartments has reached a steady state in the propagation of virus from type 1 at this time, while virus from type 2 has not yet started to spread. When type 2 virus appears, it will cause the susceptible fraction to decrease at t G [176, 215], the exposed fraction to increase at t G [180, 211], to decrease at t G [211, 240], and the infected fraction to increase at t G [208, 236], and at t G [236, 446] decrease. Dead and recovery fraction continue to rise steadily.

Fig. 11. Different time: a) only type 1 in the system; b) only type 2 in the system; c) type 1 and type 2 in the system

Table 4. Fraction of each state at the peak for different case

_E_I

Type 1 and Type 2 for same time 0.9738 0.5725/0.5679 Type 1 and Type 2 for different time 0.3099/0.6548 0.2207/0.5593

Then, we compare Figure 7c and Figure 11c, trying to find differences between the two types of virus at the same time and at different times.

Table 4 shows that when the two types of virus are present on the system at different times, the proportion of exposure and infection is smaller than when they are present on the system at the same time, and the sequential presence of the virus moderates the propagation process.

Table 5. Time required for each state to stabilize for different case

S I E D R~ Type 1 and Type 2 for same time 42 300 100 343 361 Type 1 and Type 2 for different time 215 477 250 573 553

Table 5 shows that when two kinds of virus appear in the system at different times, the time required for each state to stabilize is longer than when they appear in the system at the same time, and the successive appearance of the vires delays the overall spread of the system end time.

Next, we changed the value of one of the parameters and fixed the value of the other parameters to observe the influence of each parameter on the spread of virus. First of all, we observe the infection rate of the system when the parameter values are different. The values of are 0.3,0.5,0.7,0.9 and the values of are 0.2, 0.4,0.6, 0.8.

a) b)

Fig. 12. The infection fraction: a) corresponding to different values; b) corresponding to different u1 values

Table 6. The peak situation of the infection fraction when the value of is different

t e [0,200] t e (200,600]

M11 = 0.3 / 0.7914

ßn = 0.5 / 0.7118

ßn = 0.7 0.327 0.5855

ßii = 0.9 0.488 0.3881

Combining Figure 12 and Table 6, at t G [0,200], the larger the value of ^11, the larger the infection fraction of the system. So in order to prevent the spread of virus, we should lower the value of ^11, for example, we need to isolate susceptible compartments and exposed compartments. At t G (200,600], this pattern is no longer met because type 2 virus also participates in the propagation behavior. At this time, the value of ^22,^12 and ^21 will directly affect the change of the system infection score.

Table 7. The peak situation of the infection fraction when the value of w1 is different

t e [0,200] t e (200,800]

Ui = 0.2 / 0.7903

U1 = 0.4 0.5632 0.3339

U1 = 0.6 0.5935 0.3211

Ui = 0.8 0.6111 0.3155

Combining Figure 12 and Table 7, the parameter has a similar effect on the infection fraction as ^11, at t G [0,200], the larger the value of w11, the larger the infection fraction of the system. So in order to prevent the spread of virus, we should

lower the value of wn, that is, reduce the contact between exposed compartments and infected compartments. When w1 = 0.2, there is only one peak, Because the value of wi is too small. At t G [0,200], the infection fraction grows extremely slowly and has not yet reached its peak. The spread of type 2 virus has already begun, and the infection fraction will rise to another peak.

We then observe the exposed fraction of the system when the parameter values are different. The value of w1 is 0.2,0.4,0.6,0.8.

E-t

0 lOO 20O 300 400 SOD 600 TOO 800 Time

Fig. 13. The exposed fraction corresponding to different u1 values

Table 8. The length of time that the exposure score dropped rapidly after the first peak

length of time

Ui = 0.2 122

u 1 = 0.4 59

Ui = 0.6 39

ui = 0.8 26

Combining Figure 13 and Table 8, the larger the value of w1, the faster the drop in the later stages of type 1 virus propagation. This means that exposed compartments are converting to infected compartments at a faster rate. To hinder the spread of virus, we must minimize the rate of contact between exposed compartments and infected compartments.

We then observe the death fraction of the system when the parameter values are different. The values of are 0.3,0.5,0.7,0.9 and the values of w1 are 0.2,0.4, 0.6,0.8.

Table 9. The death fraction when the value of is different

t = 200 t = 500

ßii = 0.3 0.004031 0.8235

ßii = 0.5 0.02569 0.797

ßii = 0.7 0.2394 0.8915

ßii = 0.9 0.4102 0.9351

Table 10. The death fraction when the value of w1 is different

t = 150 wi = 0.2 0.002828 wi = 0.4 0.3865 wi = 0.6 0.4197 wi = 0.8 0.4346

Combining Figure 14 and Table 9, at t G [0,200], the larger the value of ^11, the higher the death fraction. After that, due to the spread of type 2 virus, the death fraction no longer conforms to this pattern.

Combining Figure 14 and Table 10, the larger the value of c1, the higher the death fraction.

We then observe the recovered fraction of the system when the parameter values are different. The values of c are 0.1,0.3,0.5,0.7.

Fig. 15. The recovery fraction corresponding to different w1 values

Table 11. The situation of the fastest growing stage of the recovery score

Required time Corresponding recovery fraction

UJ1 = 0.1 406 0.09131

UJi = 0.3 224 0.2077

uJi = 0.5 150 0.3973

uJi = 0.7 106 0.4896

Combining Figure 15 and Table 11, the larger the value of uj1, the faster the recovery fraction increases, and when a steady state is reached, the higher the fraction.

Among the COVID-19 outbreaks in 2019(Hannah and Ritchie, 2020), France had its first confirmed case in Bordeaux on January 24, 2020, and India had its first confirmed case on January 30, 2020. According to the data, on March 1, 2020, the cumulative number of infections per 1 million people in France was 6.32, and the cumulative number of deaths per 1 million people was 0.1. On May 1, 2020, the cumulative number of deaths per 1 million people in India was 0.1. The cumulative number of infections was 26.47, and the cumulative number of deaths per 1 million people was 0.87. On September 17, 2022, France currently has 395398 infections, 154672 deaths, 34310015 recoveries, a mortality rate of 4%, and a recovery rate of 98.4%; India currently has 46848 infections, 528302 deaths, and recoveries. was 43953374, the mortality rate was 1.2%, and the recovery rate was 98.7%.

In COVID-19, because the recovered refers to the patients who have recovered from the infection of the new type of coronavirus pneumonia. So here we assume that the rate of conversion of susceptible to cured is 0. By calculation, we can derive the average infection rate, which we assume the virus spreads at a uniform rate. Here, we only show outbreaks in France and India, where the virus emerged 60 days apart. The state transition is shown in Figure 16. In figure a), after t=200, the state change tends to be stable, in figure b), the state change has a significant change after t=600, in figure c), when 500 < t < 800, the state changes significantly .

a) b) c)

Fig. 16. a) State transition of France; b) State transition of India; c) State transition of two countries

The graph shows that the peak of virus transmission in France is July 2020 and that in India is August 2022. But when the two countries influence each other, there are no clear two peaks, mainly because the total population of France is much smaller than that of India, and the individual changes are minimal to the whole.

Liu Xiuxiu, Elena Gubar Table 12. Contrast infection fraction

Highest infection fraction Corresponding moment France 0.8779 98

India 0.5323 842

Both 0.6955 686

In Table 12, we can clearly see that the virus infection situation in the two countries is very different. When the two countries are influenced by each other due to population movement, the infection rate and the peak period of infection are between the highest and lowest values before the influence.

6. Conclusion

In this work, we have done:

1) On the epidemic base model SIR, add some new compartments, consider systems with multiple types, and model the modified SEIDR system and apply it to the spread of virus.

2) The optimal control problem is formulated in the context of a total cost minimization model. We solve this nonlinear optimal control problem by dynamically choosing the vaccine propagation rate using the Pontryagin maximum principle.

3) We focus on considering only two different types of situations in the system. We fit the propagation process of this model with Matlab. Here, we consider this an uncontrolled model. Compare two types of virus at the same time and at different times. Then change some parameters of the model and observe the effect of the parameters on the virus propagation process. Numerical computations reveal many interesting behavioral strategies.

Appendix

1. Proof of Lemma 1

Lemma 1. For any i, is a decreasing function of t.

Proof: ^i(t) is everywhere continuous (due to continuity of the states and adjoint functions) and differentiable, if ui(•) is continuous. For any t where ui(•) is continuous, we have:

U = (Af - AR)Si + (Af - AR)Si + 0i[(Xf - AR)Ii + (A' - AR)ii]. (24)

Uli

Therefore, after some regrouping and cancellation of terms, at any t, we have

= -(AE - AR)^ Si^jiEj - OiliEiUi) - (AD - AR)diIiSi Ui (25)

- (Si + °iIi)(^LrT - Yi(ui,t)) - °iIi fr.

Now, since 0 < 6i < 1, during virus propagation, have YlMjLi Si^jiEj > IiuiEi > QiIiEiui, Lemma 2 will imply that the sum of the first and second terms will be negative. The third and last terms will be non-positive due to the definitions of f (•) and L(•). So <j>i(t+) < 0 for all t. The proof for 4>i(t-) < 0 is exactly as above. In

a very similar fashion, it can be proved that <i (t) < 0 at all points of continuity of «(•), which coupled with the continuity of <i(t) shows that it is a decreasing function of time.

2. Proof of Lemma 2

Lemma 2. The positivity constraints Af (t) - Af (t) > 0 and Af (t) - Af (t) > 0 hold for all i = 1,..., M and all t G [0, T).

Proof: First, we note that Af (T) =0 at any t G [0,T] at which u is continuous, Af = — Jfr = — TTcT^ < 0. Thus, since u is piecewise continuous, Af (t) > 0 for all 0 < t < T. And by the same token, Af (t) > 0 for all 0 < t < T.

We note that Af(T) =0 at any t G [0,T] at which u is continuous, Af = ^Lff - — hi(ui) + (Af — Af )SitJiUi + (AI — Af )IiaJiUi0i = ^Lff1 — Yi(ui, t) > 0. Thus, since u is piecewise continuous, Af (t) < 0 for all 0 < t < T.

So, we can get Af (t) — Af (t) > 0 and Af (t) — Af (t) > 0 hold for all i = 1,..., M and all t G [0,T).

References

Bichara, D., Iggidr, A. and Sallet, G. (2014). Global analysis of multi-strains SIS, SIR and MSIR epidemic models. J. Appl. Math. Comput., 44, 273-292. https://doi.org/10.1007/s12190-013-0693-x Chiang, Alpha C. (1992). Elements of Dynamic Optimization. New York: McGraw-Hill Eshghi, S., Khouzani, M. H.R., Sarkar, S. and Venkatesh, S.S. (2016). Optimal Patching in Clustered Malware Epidemics. In: IEEE/ACM Transactions on Networking, 24(1), pp. 283-298. https://doi.org/10.1109/TNET.2014.2364034 Grass, D., Caulkins, J., Feichtinger, G., Tragler, G. and Behrens, D. (2008). Optimal Control of Nonlinear Processes: With Applications in Drugs, Corruption, and Terror. https://doi.org/10.1007/978-3-540-77647-5 Gubar, E., Kumacheva, S., Zhitkova, E. and Porokhnyavaya, O. (2016). Impact of Propagation Information in the Model of Tax Audit. In: Petrosyan, L., Mazalov, V. (eds.). Recent Advances in Game Theory and Applications. Static & Dynamic Game Theory: Foundations & Applications. Birkhauser, Cham. https://doi.org/10.1007/978-3-319-43838-2-5

Gubar, E., Policardo, L., Carrera, E. and Taynitskiy, V. (2021). Optimal Lockdown Policies

driven by Socioeconomic Costs. Papers 2105.08349, arXiv.org. Gubar, E. and Zhu, Q. (2013). Optimal control of influenza epidemic model with virus mutations. 2013 European Control Conference (ECC), pp. 3125-3130. https://doi.org/10.23919/ECC.2013.6669732 Khouzani, M. H.R., Sarkar, S. and Altman, E. (2011). Optimal control of epidemic evolution. 2011 Proceedings IEEE INFOCOM, 2011, pp. 1683-1691. https://doi.org/10.1109/INFCOM.2011.5934963 Khouzani, M. H.R., Sarkar, S. and Altman, E. (2012). Maximum Damage Malware Attack in Mobile Wireless Networks. In: IEEE/ACM Transactions on Networking, 20(5), 1347-1360. https://doi.org/10.1109/TNET.2012.2183642 Li, Y., Hui, P., Jin, D., Su, L. and Zeng, L. (2011). An optimal distributed malware defense system for mobile networks with heterogeneous devices. 2011 8th Annual IEEE Communications Society Conference on Sensor, Mesh and Ad Hoc Communications and Networks, pp. 314-322. https://doi.org/10.1109/SAHCN.2011.5984913 Mathieu, E., Ritchie, H., Rodes-Guirao, L., Appel, C. Giattino, C., Hasell, J., Macdonald, B., Dattani, S., Beltekian, D., Ortiz-Ospina, E. and Roser, M. (2020). Coronavirus Pandemic (COVID-19). Available at https://ourworldindata.org/coronavirus (accessed 01.11.2022).

Pastor-Satorras, R., Castellano, C., Van Mieghem, P. and Vespignani, A. (2015). Epidemic processes in complex networks. Reviews of Modern Physics, 87(3), 925. https://doi.org/10.1103/RevModPhys.87.925 Shakkottai, S. and Srikant, R. (2008). Peer to Peer Networks for Defense Against Internet Worms. IEEE Journal on Selected Areas in Communications, 25(9), 1745-1752. https://doi.org/10.1109/JSAC.2007.071212 Taynitskiy, V., Gubar, E., and Zhu, Q. (2019). Optimal Impulse Control of SIR Epidemics Over Scale-Free Networks. In: Song, J., Li, H., Coupechoux, M. (eds.). Game Theory for Networking Applications. EAI/Springer Innovations in Communication and Computing. Springer, Cham. https://doi.org/10.1007/978-3-319-93058-9-8