Decision making procedure in optimal control problem for the Sir model Текст научной статьи по специальности «Математика»

Gubar Elena and Zhitkova Ekaterina

St.Petersburg State University,

Faculty of Applied Mathematics and Control Processes,

Bibliotechnaya pl. 2, St.Petersburg, 198504, Russia E-mail: alyona.gubar@gmail.com zhitkovakaterina@mail.ru

Abstract In this work we join on classical SIR model to describe influenza epidemic in urban population with procedure of making decision. We suppose that agent in urban population makes a choice: whether or not to participate in vaccination company. Each decision involve different costs and indirectly influence on the population state. We formulated an optimal control problem to study the optimal behavior during epidemic period and vaccination company. All theoretical results are also supported by the numerical simulations.

Keywords: SIR model, vaccination problem, evolutionary games, optimal control, epidemic process.

1. Introduction

Originally Susceptible-Infected-Recovered model and its modification describe a fast spreading process, such as influenza epidemic or other forms of respiratory viral diseases, circulated in urban population. Total population is divided into three subgroups: Susceptible, Infected and Recovered. Susceptible is group where people are not infected, Infected is a group of people having the disease, and Recovered is group, where all members have immunity to the disease. Human population meets influenza epidemic almost every year then SIR model is very actual and can be used in social and economic applications.

During the years many medical methods such as preventive measures, intensive treatment, etc. were developed to protect entire population during annual epidemics. Hence preventive measures or medical treatment can be considered as an external influence to the development of epidemics and can be used as control parameters in the model. Since vaccination is one of the most effective method to protect population from the annual epidemics then we chose is as an control parameter.

However vaccination can not be absolutely effective and moreover such as it was proofed in the previous research total vaccination is very expensive and usually do not apply to protect population against a flu epidemics. Hence we can establish a new problem such as vaccination problem.

We assume that vaccination company occurs before the seasonal epidemic begins, because it is necessary take into account regulation immune system of individual after vaccination, because failing health after vaccination not allow to resist against another viruses. Unfortunately flu vaccines are effective only for one season owing to mutation of pathogens and waning immunity. We suppose that influenza epidemic continues until there are no more newly infected individuals.

In this model we assume that before epidemic period each agent of population may choose a behavior: participate or not in vaccination company. All decision provoke corresponding costs and influence on the future agent’s incomes. If agent of population chooses the application of vaccination then he should pay it costs and estimate the consequences if vaccine is not effective. Hence in current work we extend classical Susceptible-Infected-Recovered (SIR) model with the procedure of making decisions. We reformulate original model in terms of optimal control and couple it with the process of making decisions.

1.1. Related Works

Recent literature has seen a large amount of interest in using optimal control and game-theoretic methods to study disease control of influenza for public health. First, this research problem was refereed in (Kermack and Kendrick, 1927), where an Susceptible-Infected-Recovered model has been proposed to study the epidemic spread in a homogeneous population. It provides a deterministic dynamical system model as the mean field approximation of the underlying stochastic evolution of the host subpopulations. In (Behncke, 2000) and (Kolesin and Zhitkova, 2004), many variants of optimal control models of SIR-epidemics are investigated for the application of medical vaccination and health promotion campaigns. In the paper (Fu et al., 2010) the vaccination problem is considered from the individual agents’ point of view.

Also epidemic models can be applied to the different fields of human activity, for instance in (Khouzani et al., 2010; Khouzani et al., 2011), optimal control methods have been used to study the class of epidemic models in mobile wireless networks, and Pontryagin’s maximum principle is used to quantify the damage that the malware can inflict on the network by deploying optimum decision rules.

Different from the work done in the past, in current work population agents based on the available information make a decision whether or not they participate in vaccination company. Their choices are included in the SIR model and as a result we receive the optimal control strategy (intensity of vaccination) that depends on the decision procedure.

2. Model

We use Susceptible-Infected-Recovered model to describe epidemiological process in urban population with following assumption that each agent in population allows to participate in vaccination company or refuses it. In current model vaccination company establish the influence to the population and hence we can consider it as control parameter in the model. Then, at time t, ns, nj, nR correspond to fractions of the population who are susceptible, infected and for all t, condition N = ns +nj+nR is justified. Define

S(t) = ^r, m = R(t) = (R(t) = 1 - S(t) - I(t j)

as portions of the susceptible, the infected and the recovered in the population.

And in addition to the above the model is formulated as follows (Khatri, 2003; Kermack and Kendrick, 1927):

dS T

—— = —dbl — u;

£ (1)

— = SSI — al; dt ’

here transmission rate from state S to I is

6 = 60m(^) = 60ml, (2)

where value S0 is a transmissibility of disease, m is a number of contacts per time unit, and parameter a = ^ which is intensive rate of transition from infected to recovered and it is in inverse proportion to the average duration of the disease. Here R = 1 — S — I, variable u(t) G (0,1) is control parameter which is interpreted as the intensity of vaccination in agents per day.

2.1. Objective Function

In this work we will minimize aggregated cost in time interval [0,T], hence at any given t following costs exist in the system: fi(I(t)) these are individual’s treatment costs, which are non-decreasing and twice-differentiable, convex functions, such as /i(0) = 0, fi(I(t)) >0, i = 1, N for I(t) > 0; functions are agent’s benefit

rate, which arise when infected agent becomes recovered, li(R(t)) is non-decreasing and differentiable function and l(0) = 0; functions hi(u(t)) describe vaccination costs that help to reduce epidemic spreading, hi(u(t)) is twice-differentiable and

increasing function in Ui(t) such as /i*(0) = 0, hi(x) >0, i = l,N when u > 0.

Hence costs function for i-th agent in population is:

Ji = fi(I (t)) — li(R(t)) + hi(u(t)). (3)

Therefore aggregated system costs is:

T N

J = \ J2(fi(I(t)) — li(R(t)) + hi(u(t)))dt. (4)

0 i=1

2.2. Making decision procedure

In current section we present a procedure of making decisions that influence to the epidemic process in urban population. Previous researches have proofed that vaccination company as a preventive measure is very effective and allows to reduce the quantity of infected in entire population. However each agent in population have a possibilities to estimate his own profit of participation in vaccination company. Agent can take into account the vaccination cost, feasible complications after vaccination and also he can estimate the herd immunity. We suppose that the last circumstance does not presume than an agent necessarily knows the exact information, he can evaluate the average number of his contacts, the current epidemic situation, that can be presented in mass communication media, etc. Meanwhile the collective result of vaccination decisions determines the level of population immunity and the strain of the epidemic in current period. When level of vaccination coverage in total population is increased then even agents who are unvaccinated have less risk to become infected. Then we assume that every agent, having this information might decide to decline the vaccination this year and thereby he reduces own vaccination costs. However agents might have incomplete information

with the some rumors from the neighbors or friends, or they also may estimate the epidemic situation incorrectly, thus this scenario leads following problem, increasing of unvaccinated individuals provoke the diminution of the herd immunity in the future and thereby collective costs during epidemic period will arise. The reduction of the vaccinated individuals induces the increasing of infected in population the then it leads that the frequency of meeting with infected agents is also increased. Then each unvaccinated agent may transform to the infected and then he should pay treatment costs, that include healthcare expenses, lost productivity and the possibility of pain. Usually treatment costs exceed the vaccination expanses.

Thus in current work we suppose that each agent chooses between two possible alternatives:

— to be vaccinated;

— not to be vaccinated and probably to be infected;

If agent participate in vaccination company then he gets a vaccination costs, in our model these costs are described by functions hi(u), where u is intensity of vaccination. Vaccination costs contain the immediate monetary cost, the opportunity cost of time spent to get the vaccine and any health effects. We also suppose that vaccination is not absolutely effective and vaccination company should be finished before the epidemic starts.

Infected agents incur treatment costs, which are denoted as functions fi(I), and when agent convalesce then his treatment costs are reduced to the value l(R), which is benefit function.

Then describe the decision procedure, each season an agent adopts one of the alternative, which determines whether or not he vaccinated. At the end of the season each agent decides whether to change the vaccination decision or not, depending on the current aggregated costs. Then agent i selects at random agent j, and in imitates his role model if opponents payoff is higher. Define probability that agent i adopts behavior of agent j as follows (Fu et al., 2010):

where pj is agent’s payoff on j-th decision, parameter j3 G (0, to).

We incorporate this probability to the basic Susceptible-Infected-Recovered model, which is presented in section 2., thus transmission rate from S to I can be rewritten:

3. Structure of optimal control

We use Pontryagin’s maximum principle (Pontryagin et al., 1962), to find the optimal control u = (ui,u2) to the problem described above in Section 2.. Define the associated Hamiltonian H and adjoint functions As, A^, AIr, AR as follows:

1

(5)

Pt3 1+exp {-l3{Pj-pi)Y

S = SomIpij.

(6)

N

H = — A0 ^2 (fi(I(t)) — li(R(t)) + hi(u(t))) +

As (-SS(t)I (t) — u) + AT (SS(t)I (t) — aI(t))

N

— A0 (fi(I(t)) — li(R(t)) + hi(u(t))) —

SS(t)I (t)(As — AT) — As u — AT a.

We construct adjoint system as follows:

N

Xi(t) = = A0 Xj + 5S{Xs{t) - Ai(t)) + Xi(t)a; ^

*s(t) = -|f = As\t)SI(t) - X2(t)SI(t) = SI(t)(Xs(t) - A/(t)); with the transversality conditions given by

Ai(T) = 0, As(T) = 0, Ar(T) = 0. (9)

According to Pontryagin’s maximum principle, there exist continuous and piecewise continuously differentiable co-state functions Ai that at every point t G [0, T]

where ui and u2 is continuous, satisfy (8) and (9). In addition, we have A(t) =

(A0(t),As(t),AI(t),AR(t)) oaeafi reoi

uGarg max H(X, (S, I, R),u). (10)

«£[0,1]

To determine an optimal control parameter that maximize Hamiltonian (7) we

, • • dH

consider derivative ——: du

N N

= _^° $3 _ As = -(A0 ^2 K(u) + As)- (11)

i=1 i=1

Now let equal to zero right parts of equations (11), (8):

N

—(A0 J2 hi(u) + As) = °; i=1

N

A0Z f(I (t)) + SS(As (t) — Ai (t)) + Ai (t)a = 0;

SI(t=1(As (t) — Ai (t))=0.

From the first equation of system (12), Hamiltonian reaches maximum if and only if next condition is satisfied:

(12)

N

1 'Su)

i=1

Let be A0 = 1, then expression (13) can be reformulated:

(A0J2 h'i(u) + As) < 0. (13)

N

Y^hi(u) < —As, (14)

1

and we will proof that As < 0. From (12) we received that

1 n

as =—£//(/(*)), (15)

rr < *

a i=1

where a > 0, ^ f'(I(t)) > 0 by definition then As < 0, hence maximum i=1 i

of Hamiltonian is reached on the negative half-space then we should proof that

function As is increasing. Consider adjoint system (8) and show that derivative

As = SI(t)(As(t) — Ai(t)) > 0.

We will proof this statement base on the next two properties (Khouzani et al., 2011):

Property 1. Let w(t) be a continuous and piecewise differential function of t. Let w(t1) = L and w(t) > L for all t G (t1,... ,t0]. Then w(t+) > 0, where w(t+) = lim v(x).

X —— Xo

Property 2. For any convex and differentiable function y(x), which is 0 at x = 0,

y'(x)x — y(x) > 0 for all x > 0.

Step I. Consider instant time moment t = T, from transversality conditions

n

(9) we have As(T) — Ai(T) = 0, and As(T) — Ai(T)= — £ f'(I(T)) < 0, Ai(T) =

i

n

Y,f'(I(T)) > 0, therefore function AI is increasing on the interval [0,T].

i

Step 2.(Proof by contradiction).

Let 0 < t* < T be the last instant moment at which one of the inequality constraints are performed:

Condition 1. AI (t) > 0, As (t) — AI (t) = 0 for t* <t < T.

Condition 2. AI (t) = 0, As (t) — AI (t) < 0 for t* <t < T.

Now consider a difference:

N

As (t*+) — Ai (t+*) = SI(T )(As — Ai ) — (A0 £ f!(I (t)) + SS(As — Ai ) + Ai a) =

i=1

6I(T)(XS - A/) - X0J2 fi(I) + T + ~fu + ~TaI ~ A/Cr i=1 11 1

AN

+ f(E(/*(/(t)) - M-^)) + hi{u(t)))) =

i=1

A N N H

= si(T)(xs - a,) mm - E №))) + j I i=1 i=1 I

E ^(-R(t)) + E + -j~u + YaI ~ A/Cr

i=1 i=1

(16)

The system ODE is autonomous, i.e., the Hamiltonian and the constraints on the control u do not have an explicit dependency on the independent variable t. Then at time t = T Hamiltonian is:

N

H (T) = —A0E(fi(I(T)) + li(R(T)) + hi(u(T)))— (17)

—SlS(T )I (T )(As (T) — Ai (T)) — As (T )u(T) — Ai (T )aI(T).

costs functions follow the next conditions fi(I(T)) > 0, li(R(T)) > 0, hi(u(T)) > 0 and transversality conditions (9) at time moment T are justified then

H(T) < 0.

Hence as far as functions f, l, h are non-decreasing we have:

N

H(t) - Ao £ (fi(I(t)) + li(R(t)) + hi(u(t))) = (1g)

-5S(T)I(T)(As(T) - A,(T)) - As(T)u(T) - A,(T)aI(T) < 0.

A N N

By property 2. the term is nonnegative -j- E wWV -£ fi(1 w > 0,

i=l i=l

from condition 1 (AS - A,) = 0 and A, > 0 and from (18) we received that As(t*+) ~ Ai(t+*) < 0, then ^(As(t*+) ~ A/(t*+)) < 0, which contradicts property 1, thus time moment t*+ does not exist.

Step II. Consider formula (16) and suppose that condition 2. is satisfied, by

A N N

property 2. we have -y- E m o))1 -E fi(I(t))) > 0, for all t it is justified that

i=1 i=1

H(t) < 0, therefore ^(\s(t*+) - A/(t*+)) < 0. This also contradicts property 1. and then time moment t*+ does not exist in case II. Hence for all t G [0,T] condition ^(As(t) — Ai(t)) > 0 is satisfied.

These proofed results can be formulated as lemmas.

Lemma 1. For all t, 0 < t < T following conditions hold (AS(t) - Ai(t)) > 0 and A, (t) < 0.

Lemma 2. For all t on time interval 0 <t <T we have

N

(.AoE hi(u) + As) < 0. (19)

i=1

Based on previous research (Khouzani et al., 2010; Khouzani et al., 2011; Pontryagin et al., 1962), we show that an optimal control u(t) = (u1(t),u2(t)) has following form.

Theorem 1. Optimal control program u(t) has following structure:

For all t such as 0 <t <t*, u(t) satisfies:

A N N

(2°)

° i=1 i=1

For all t such as t* < t < T:

u(t) = 0.

4. Numerical simulations

In this section we present numerical simulation which are used to illustrated the structure of the optimal control and influence of the human decision to the epidemic process. In the example we suppose that population size is N = 1000, initial fraction of subpopulations are: S(0) = 950, I(0) = 50, R(0) = 0. Following values of the system parameters are used in the simulation: h = 0,1 is model step, £o = 0,06 is

transmissibility of decease, l = 35 is number of contacts per time unit, a = 1/15 = 0,06(6) is intensity of recovering.

Based on previous research and on the result of the experiments in this work we assume that typical duration of disease is 10-15 days. After simulation we receive that with mentioned initial states and auxiliary parameters the maximum quantity of infected is I(t) = 511 and epidemic peak is reached at t = 15 day.

In figure 1 Susceptible-Infected-Recovered model is presented:

Fig. 1: SIR model without application of the control

Below we present the case, where we apply optimal intensity of vaccination (optimal control strategy), which allow to reduce the number of infected in population. In current model vaccination was used as a control parameters in the system hence agents from subpopulation Susceptible directly transfer to subpopulation Recovered, obtaining immunity. After numerical simulations for the same initial data we get that the maximum quantity of infected is I*(t) = 426 at time t = 13. Therefore we can see that maximum number of Infected is less than in previous case and in comparing with Fig.1 epidemic peak is achieved early.

Fig. 2: Application optimal control to the SIR model.

In figure 3 we illustrate the optimal control that minimize aggregated costs on the preventive measures. For the considered initial states optimal control will be switched off at t* =7 day.

One of the main aim of the work is to show that participation of agents in vaccination company reduces aggregated costs of the entire population, hence in following figures we present aggregated costs for different cases. First individual

u

10 20 30 40 50

Fig. 3: Optimal control, uopt = 15 agents per day.

agents costs are define as follows: treatment costs are fi(I(t)) = aI(t) + b, where a =1, b = 1.4, vaccination costs hi(u(t)) = k1u2 + k2, and k1 = k2 = 0.5, li(R(t)) = cR(t) + d, c =1, d = 0.05.

In figures 4.-5 aggregated costs received for the time interval [0, T] and they are equal to J(u) = 801.5 monetary units (m.u.), which can be in US dollars, Rubles or Euros depending on the context.

Aggregated costs

Fig. 4: Aggregated system cost for SIR model without application of the optimal control.

If optimal control is applied to the system then aggregated system costs decrease and value of functional is equal to J(u) = 707,14 m.u.

Fig. 5: Aggregated cost for SIR model with application of the optimal control.

To complete the our illustrative example let’s consider a modification of the model, where we take into account only agents choices. Let’s say that for instance

that intensity of vaccination is u = 5 agents per day, it means that only five agent accept a decision about the participation in vaccination company. We should add also that u < uopt. In such case dynamics in SIR model is changed, the maximum number of infected is achieved at t* = 14 day and equal to Imax = 482 agents. The result of simulation is presented in figure 6.

I

0.»

o.« \ ,

Ok» \ /*

0.» N / \ ✓’

0,5

0.4 / X /

0J / V* / * \ \

0.2 1 / . N \ S * ^

04 1/ J,* v. —

0 _ m.j |20 n 10 40 50

Fig. 6: SIR model with control parameter and decision making procedure.

We can see that the quantity of infected exceeds the number of infected in situation, where we apply optimal control strategy to the population, and the number of agents, choosing vaccination decision is not enough to protect population during the epidemic season. Hence in this case we also can show that aggregates costs increase and value of functional is J(u) = 761.36 m.u.

Aggregated costs

2,5

2 -

1,5 -

1 -0,5 -

10 20 30 40 50

Fig. 7: Aggregated costs with u = 0.005 is equal to J(u) = 761.36 m.u.

To summarize the results we have to add that since agents in population do not have reliable information about epidemic situation then they may choose incorrect decision, which provoke the degradation of epidemic state in total population.

5. Conclusion

In this paper, we have studied an epidemic model that takes into account the agent motivation to participation in the vaccination company. We incorporate procedure of making decision to the simple Susceptible-Infected-Recovered model and have formulated this model in special case. Using Pontryagin’s maximum principle, we have shown the structure of optimal control, which is depending on the agents costs induced by choosing decisions. We supported our results with numerical simulations,

observing different cases of epidemic process in entire urban population. In future work we would extend this model including different structure of population, it means that human decision may depend on his social group, not only his costs and to modify the model, using number of contacts as a function of the time.

References

Behncke, H. (2000). Optimal control of deterministic epidemics. Optim. Control Appl. Meth., 21, 269-285.

Bonhoeffer, S., May, R. M., Shaw, G. M. and M. A. Nowak (1997). Virus dynamics and drug therapy. Proc. Natl. Acad. Sci. USA, 94, 6971-6976.

Conn, M. (2006). Handbook of Models for Human Aging. Elsevier Academic Press, London (ed.).

Fu, F., D. I. Rosenbloom, L. Wang and M. A .Nowak (2010). Imitation dynamics of vaccination behaviour on social networks. Proceedings of the Royal Society. Proc. R.Soc. B., 278, 42-49.

Gjorgjieva, J, Smith, K., Chowell, G., Sanchez, F., Snyder, J., and C. Castillo-Chavez (2005). The role of vaccination in the control of SARS. Mathematical Biosciences and Engineering, 2(4), 753-769.

Gubar, E., Zhitkova, E., Fotina, L., and I. Nikitina (2012). Two Models of the Influenza Epidemic. Contributions to game theory and management, 5, 107-120.

Karpuhin, G. I. (1986). Gripp . Leningrad.: Medicina (ed).

Khatri, S., Rael, R., J. Hyman (2003). The Role of Network Topology on the Initial Growth Rate of Influenza Epidemic. Technical report BU-1643-M.

Kermack, W. O. and A. G. Mc Kendrick (1927). A contribution to themathematical theory of epidemics. In: Proceedings of the Royal Society. Ser. A., Vol. 115, No. A771, pp. 700-721.

Khouzani, M.H.R., Sarkar S. and E. Altman (2011). Optimal Control of Epidemic Evolution. In: Proceedings of IEEE INFOCOM.

Khouzani, M. H.R., Sarkar, S. and E. Altman (2010). Dispatch then stop: Optimal dissemination of security patches in mobile wireless networks. Proc. of 48th IEEE Conference on Decisions and Control (CDC), pp. 2354-2359.

Mehlhorn, H. et al. (2008). Encyclopedia of Parasitology. Third Edition Springer-Verlag Berlin Heidelberg, New York (eds).

Pontryagin, L. S., V. G. Boltyanskii, R. V. Gamkrelidze et E. F. Mishchenko (1962). The Mathematical Theory of Optimal Processes. Interscience.

Sandholm, W. H., E. Dokumaci and F. Franchetti (2010). Dynamo: Diagrams for Evolutionary Game Dynamics, version 0.2.5. http://www.ssc.wisc.edu/ whs/dynamo. Kolesin, I.D., E. M. Zhitkova (2004). Matematicheskie modeli epidemiy. SPbU.