Определение точности математического моделирования характеристик эпидемии гриппа

13. Ярошенко Е.Б., Бурневич Э.З., Мойсюк Я.Г. Роль вирусных гепатитов в развитии гепатоцеллюлярной карциномы // Практическая онкология. 2008. № 4 (9). С. 189 - 193.

14. Готье С.В., Мойсюк Я.Г., Ибрагимова О.С. Органное донорство и трансплантация в Российской Федерации в 2006 - 2010 годах: Третье сообщение регистра Российского трансплантологического общества // Вестник трансплантологии и искусственных органов. 2011. Т. XIII. № 2. С. 6 - 20.

15. Готье С.В., Мойсюк Я.Г., Ибрагимова О.С. Органное донорство и трансплантация в Российской Федерации в 2009 году: Второе сообщение регистра Российского трансплантологического общества // Вестник трансплантологии и искусственных органов. 2010. Т. XII. № 3. С. 6 - 15.

16. Материалы II Конференции Межрегиональной общественной организации «Общество трансплантологов» (тезисы докладов) // Трансплантология. 2010. № 1. С. 72 - 114.

17. Протокол диагностики и лечения больных вирусными гепатитами В и С // Российский журнал гастроэнтерологии, гепатологии и колопроктологии, 2010. № 6. С. 4 - 60.

18. Козлова А.В. Профилактика и лечение возвратных инфекций после ортотопической аллотрансплантации печени у больных циррозом печени HBV- и HCV-этиологии: Автореф. дис. ... д. м. н. - М., 2011.

19. Хубутия М.Ш., Никулина В.П., Годков М.А. и др. Клинико-лабораторные аспекты монокомпонентной иммуносупрессии при трансплантации печени // Трансплантология. 2009. № 2. С. 25 - 30.

20. Полный вариант статьи будет опубликован во № 2(69) 2013 журнала «Эпидемиология и Вакцинопрофилактика».

21. Консолидированный бюджет Российской Федерации и бюджетов государственных внебюджетных фондов в 2009 году, за исключением акцизов по подакцизным товарам (продукции), доходов от использования имущества, находящегося в государственной и муниципальной собственности, безвозмездных поступлений (http://www.gks.ru/bgd/regl/b10_11/IssWWW.exe/Stg/d2/23-01.htm).

22. http://www.consultant.ru/online/base/?req=doc;base=law;n=93318/

Определение точности математического моделирования характеристик эпидемии гриппа

Б.М. Десятков

ФБУН «ГНЦ ВБ «Вектор» Роспотребнадзора, пос. Кольцово Новосибирской области

Резюме

В работе исследуется точность теоретических результатов, полученных с использованием математической модели, путем сравнения их с данными по суммарной заболеваемости населения в 45 городах России во время эпидемии гриппа A(H1N1)v и сезонной эпидемии гриппа и ОРВИ в период с сентября 2009 по январь 2010 года.

Выявлена достоверная корреляционная связь между значениями численности населения города S0 и параметром v, характеризующим вероятность инфицирования. Однако использование корреляционной связи при вычислении основных характеристик эпидемий гриппа не имело успеха. Причиной этого является большая чувствительность математической модели к неизбежным ошибкам в задании значений входных параметров модели. Для преодоления указанных недостатков математической модели предложен новый метод, в котором задаются значения не отдельных параметров v и S0, а функция f(t) = A - B • t, которая является аппроксимацией произведения v • S(t), вычисленного по данным наблюдений. Тестовые расчеты эпидемий гриппа с использованием этого метода дали удовлетворительные результаты.

Ключевые слова: математическая модель, эпидемия гриппа, точность результатов расчетов

The Evaluation of the Accuracy of Mathematical Modeling Characteristics of an Influenza Epidemic

B.M. Desyatkov

«Vector» State Research Center of Virology and Biotechnology, Novosibirsk Region, Koltsovo

Abstract

The accuracy of the theoretical results obtained using the mathematical model by comparing them with data for total morbidity in 45 Russian cities during the epidemic of influenza A(H1N1)v and the seasonal epidemic and acute respiratory viral infection in the period from September 2009 to January 2010 is investigated.

A significant correlation was revealed between the number of the city inhabitants S0 and the parameter v, which characterizes the probability of infection. However, the attempts to use this correlation for making a mathematical prediction of influenza epidemics were unsuccessful. This is due to the great sensitivity of the mathematical model to inevitable errors in setting the values of the model input parameters. To overcome these shortcomings of the mathematical model, a new method was proposed in which the function f(t) = A - B • t that is an approximation of the product v • S(t) and calculated from the observed data. Is set instead of the values of individual parameters v and S0. The test calculations of influenza epidemics using this method gave satisfactory results. Key words: mathematical model, influenza epidemic, accuracy of the calculations

Введение

Большинство современных математических моделей, описывающих развитие эпидемических процессов, являются детерминистическими, когда оце-

ниваются только средние в каждый момент времени значения ожидаемого числа инфицируемых индивидуумов и выбывших из числа восприимчивых. Однако развитие эпидемии - по существу случайный

процесс, описание которого с привлечением только средних показателей, без информации о статистическом разбросе, не позволяет объективно оценить минимальные и максимальные ресурсы, требующиеся для проведения лечебно-профилактических мероприятий.

Вследствие этого разработка статистических моделей эпидемических процессов представляется актуальной для решения задач практической эпидемиологии. Наличие таких моделей позволяет создавать различные сценарии эпидемического процесса с оценкой возможного диапазона экономического ущерба и ресурсов, необходимых для ликвидации последствий эпидемии.

В работах [2, 6, 7] описываются стохастические модели, основанные на использовании функции плотности вероятности для описания перехода от одного случайного состояния к другому. Моделируется одно из возможных сценариев развития эпидемии. Для последующего вычисления математического ожидания и дисперсии необходимо иметь большой набор таких сценариев.

В работах [3 - 5] предложена статистическая модель, позволяющая проводить моделирование эпидемических процессов с оценкой их статистических характеристик путем использования дифференциальных уравнений для дисперсии. Опыт работы с такой моделью показал, что необходима тщательная проверка исходной модели, включающая использование реальной информации о заболеваемости во время эпидемий в разных городах, с последующей ее доработкой.

Цель данной работы - проверка точности теоретических результатов, полученных с использованием математической модели, путем сравнения их с данными по суммарной заболеваемости населения в 45 городах России во время эпидемии гриппа A(H1N1)v и сезонной эпидемии гриппа и ОРВИ с сентября 2009 по январь 2010 года.

Материалы и методы

Работа выполнялась в два этапа. На первом этапе анализировались только реальные данные по недельной заболеваемости гриппом с целью обнаружения корреляционных связей между основными характеристиками эпидемий и использования их при моделировании эпидемического процесса. Данные по заболеваемости получены в Роспотребнадзоре. Представлены: зависимость максимального числа заболевших P от числен-

max

ности населения городов S0; зависимость времени появления максимального числа заболевших T(Pmax) от численности населения городов S0; зависимость временн го интервала c момента быстрого развития эпидемии до появления максимального числа заболевших от численности населения городов S0. Изучалась связь между превышением эпидемического порога и развитием эпидемии. На втором этапе анализировались данные по заболеваемости совместно с результатами матема-

тического моделирования эпидемического процесса. Для этого использовалась известная математическая модель SIR («Susceptible» - здоровые люди, восприимчивые к заболеванию; «Infectious» - инфицированные больные; «Recovered» - переболевшие моделируемым заболеванием люди, больше к нему не восприимчивые) [2, 6, 7].

Исходная система дифференциальных уравнений решалась конечно-разностным методом. Для этого она аппроксимировалась следующими разностными уравнениями:

S

I

R(t + At) = R(t) + At • w • I(t),

где S - число чувствительных; I - число инфицированных; R - число выбывших из группы инфицированных; V - параметр, характеризующий вероятность инфицирования;

w - частота выбывания индивидуума из группы I; величина в = 1/w - длительность протекания болезни, для гриппа это 5 - 7 суток; t - время, At - шаг интегрирования по времени.

Второе слагаемое в правой части первого уравнения P = At •v• S(t) • I (t) описывает число инфицированных за временной интервал At. Поскольку мы располагаем только недельными значениями числа заболевших, то именно эта величина за недельный интервал вычислялась и анализировалась. При этом мы предполагали, что все инфицированные заболевают.

В начальный момент времени задавались значения Io- число инфицированных в первую неделю эпидемии и So - численность населения выбранного города.

С позиции математической модели особенности развития эпидемии зависят от четырех параметров: от численности населения S0, от начального числа инфицированных Io, от значений w и v. Последний параметр наиболее важен и, к сожалению, наименее изучен. Значение этого параметра определяют такие факторы, как климатические и погодные условия, экологическая ситуация, иммунный статус населения, масштаб вакцинирования, интенсивность перемещения населения и др. Роль каждого из указанных факторов на качественном уровне очевидна, однако дать им количественную оценку практически невозможно. Поэтому приходится решать обратную задачу - оценивать или подбирать значение параметра путем сравнения полученных на модели теоретических результатов с реальными наблюдавшимися (зарегистрированными) значениями основных характеристик эпидемии.

Таким образом, с использованием модели (1) были выполнены следующие исследования. Для всех городов решена обратная задача - проведено моделирование эпидемий и вычислены значения

параметров модели v и w, при которых достигнуто наилучшее совпадение рассчитанных и зарегистрированных значений характеристик эпидемии. Был рассмотрен достаточно широкий диапазон разумных значений этих параметров, но только в области указанных значений обнаружилось наилучшее совпадение. Этот факт доказывает единственность полученного решения обратной задачи.

Исследовано влияние параметров v и w на характеристики эпидемии. Построена зависимость параметра от численности населения S0. Построена зависимость параметра от максимального числа заболевших Pmax = Max(M-V • S(t) • I (t)). Исследовано влияние параметров S0 и I0 на характеристики эпидемии.

Перед математическим моделированием эпидемического процесса и сравнением полученных теоретических результатов с данными по заболеваемости населения были проведены расчеты по оценке чувствительности используемой математической модели к неизбежным ошибкам в задании входных данных. Для этого каждый из четырех параметров модели последовательно задавался с некоторой ошибкой, при этом остальные параметры оставались неизменными. Моделировалось развитие эпидемии и сравнивались максимальные значения числа заболевших за неделю в условно «точном» решении, полученном исходя из параметров, полученных при решении обратной задачи, и в «приближенном» решении. Результаты этих расчетов представлены в таблице 2.

X - X

ÖX = abs(-) -100%,

X

P -P

AP(£X) = abs(-max-max) -100%,

P

max

где:

6X - относительная ошибка при задании параметровX = I0, S0, v и w; ~ и Pmax - «приближенные» значения; AP(X) - абсолютное значение относительной ошибки при вычислении P .

1 max

Результаты и обсуждение

Результаты первого этапа работы

Типичное развитие эпидемии гриппа представлено на рисунке 1, где линия с квадратными точками - изменение со временем недельного числа заболевших в Новосибирске, линия с круглыми точками - изменение со временем показателя превышения эпидемического порога, сплошная линия - рассчитанное по модели число заболевших. Как и для большинства эпидемий, четко прослеживается связь между резким увеличением числа заболевших и превышением эпидемического порога. Как только этот показатель становится больше нуля, начинается быстрое развитие эпидемии.

Такая зависимость проявляется почти для всех городов, за исключением шести.

Наблюдается убедительная линейная связь между максимальным числом заболевших P и на-

max

чальным значением числа чувствительных S0. Соответствующая аппроксимирующая функция имеет вид: Pmx = 1,648 + 0,016 • S0 с уровнем достоверности R2 = 0,81. Но нет корреляционной связи между начальным числом инфицированных I0 и максимальным числом заболевших P , а также време-

max

нем наступления пика эпидемии. Для городов с примерно одинаковой численностью населения эти характеристики существенно различаются. Пик эпидемии наблюдается в интервале 9 - 13 недель и не зависит от численности населения S0. Величина временного интервала от момента быстрого нарастания числа заболевших до максимального значения находится в диапазоне 1 - 3 недели и тоже не зависит от S0.

Результаты второго этапа работы

В результате решения обратной задачи для каждого города получены значения параметров w и v, при которых наблюдается наилучшее совпадение вычисленных и зарегистрированных значений основных характеристик эпидемии. При этом ошибка в вычислении P не превышала 0,5%, ошибка в

max

определении момента пика эпидемии была меньше 1 недели. Был рассмотрен достаточно широкий диапазон разумных значений этих параметров, но только в области указанных значений обнаружилось наилучшее совпадение. Это доказывает единственность решения обратной задачи. В таблице 1 представлена часть этих результатов, относящихся только к первой группе из шести. Результаты для остальных пяти групп - аналогичные. Подробнее о группах будет сказано ниже. Для большинства городов продолжительность болезни - значение параметра в = 1/w - находится в интервале 4,0 - 6,5 дня. Это согласуется с общепринятой продолжительностью заболевания гриппом. Однако для девяти городов в < 3,8 дня, а для трех городов (шестая группа) он равен 2,8, 2,6 и 2,0 дня.

На рисунке 2 приведена зависимость значения параметра v от S0 (точки) - сплошная кривая, соответствующая аппроксимации с уровнем достоверности R2 = 0,93. Наблюдается очень четкая монотонно убывающая корреляция - чем больше численность населения города, тем меньше значение параметра v. Именно эта корреляционная связь была исследована далее при моделировании эпидемии. Аналогичная зависимость наблюдается между параметрами v и Pmax.

В таблице 2 приведены осредненные по всем городам значения ошибок. Из этой таблицы видно, что математическая модель эпидемического процесса оказалась менее всего чувствительна к ошибкам в задании параметра I0. Даже ошибка в 50% приводит к ошибке вычисления P не более

max

15%. Однако ошибки при задании других параме-

Рисунок 1.

Число заболевших за неделю (квадратные точки), показатель превышения эпидемического порога (круглые точки) и рассчитанное число заболевших (сплошная кривая) для Новосибирска

3

ш

ф

^

о ю со со о

о

25000 _

20 000

15 000

10 000

5000

н—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I—ь

3 5 7 9 11 13 15 17 19 Номер недели

160 140 120 4- 100 80 60 40 20 0 -20 -40

о р

о

с

д и

с

т ее

и н

е

3

-О

ш

е р

1= л

л

е

£

со со

V

о

С

Таблица 1.

Результаты анализа эпидемий в городах первой группы

Первая группа

Города ^ чел ^ чел Ртах, чел. V -10-6 (чел./день)-1 п,неделя в = 1М день А, (день)-1 В, (день)-2

Киров 457 400 1895 8977 0,3881 13 6,7 0,1851 0,0027

Кемерово 485 000 1608 13 757 0,3952 13 6,5 0,2042 0,0038

Астрахань 506 400 2014 9784 0,3699 13 6,3 0,1951 0,0028

Нижний Новгород 1 311 200 4974 23 427 0,1418 13 6,3 0,1931 0,0026

Смоленск 325 500 1306 5822 0,5957 12 6,0 0,2011 0,0027

Тверь 409 400 1296 7671 0,4861 13 5,9 0,2073 0,0028

Вологда 292 800 1251 7441 0,7087 11 5,8 0,2182 0,0040

Ульяновск 635 600 3861 11 220 0,3089 10 5,8 0,2015 0,0029

Среднее по группе 6,16 0,2007 0,0030

тров гораздо сильнее влияют на точность результата. Так, даже относительно небольшая ошибка 6у= 10%, а тем более ошибка 6у = 20% приводят к практически неприемлемым результатам.

Аналогичные выводы справедливы для оценки влияния ошибок при задании параметров 50 и V. К сожалению, именно параметр V является наименее изученным и его практически невозможно

определить с необходимой точностью. Произвести количественную оценку влияния всех этих факторов можно только при статистическом анализе большого количества эпидемий для каждого города. Параметр 50 также оценивается с недостаточной точностью. Формально это число чувствительных индивидуумов, которые потенциально могут быть инфицированы. В общем случае его значение

0

Рисунок 2.

Зависимость параметра V от численности населения S0

Точки - результат анализа, сплошная линия - аппроксимирующая функция

1-г

4000 6000

Sq, тыс. чел.

Рисунок 3.

Графики средних значений трендов Щ) = А - В  t для шести групп: нижний график соответствует первой группе, верхний - шестой группе

Таблица 2.

Абсолютные значения усредненной по городам относительной ошибки АР вычисления величины Рт при задании значений параметров S0, V и w с относительной ошибкой б

Параметры a, %

1 2 5 10 20 50

A P(6 I,) 0,3 0,6 1,6 3,3 6,0 15,1

P(S 13,1 27,2 76,5 174,8 403,0 1050,0

P(6v) 12,2 25,2 69,5 157,8 328,8 674,7

P(6w) 11,2 23,2 64,5 149,1 371,1 1256,9

Рисунок 4.

Результаты точности вычисления характеристик эпидемии гриппа; АРтах - относительная ошибка в вычислении максимального числа заболевших (величины пика эпидемии), Т (Рпаа) - ошибка в определении момента наступления пика эпидемии

80

60

40

£

J <1

-20

-40

-60

20

0

I

• j • • : -1 . 1 1

i 1 -3 -2 -• • - • « i i 1 0 ' - ' i '  • 1 • * 1 2 3 AT (Pmax), неделя

может отличаться от численности населения для конкретного города. Однако мы вынуждены задавать его равным численности населения, которое тоже точно не известно на момент возникновения эпидемии. Информация о параметре w более определенная. Величина в = 1является длительностью инфекционного периода, и, как следует из таблицы 2, ошибка в задании этого параметра в первый день (примерно 20%) также приводит к очень большим ошибкам в вычислении параметра Ртах. При задании «ложных» значений одновременно для нескольких параметров величина ошибки результата может быть еще больше. Например, при

задании параметров v и S0 с относительной ошибкой 1, 5 и 10% величина относительной ошибки AP получилась равной соответственно 16,6, 96,7

max

и 227,4%.

Очевидно, что в таких условиях трудно получить более или менее точное описание развития эпидемии, задавая значение входных параметров I0, S0, v и w. Это подтверждается результатами следующих расчетов. Предварительно по дискретным данным зависимости параметра v от численности населения S0 была построена аппроксимирующая ее непрерывная функция v = 0,2269• S-0'9935 (сплошная кривая на рис. 2) с уровнем достоверности

R2 = 0,93. Далее выполнялся расчет - выбирался город, задавались соответствующие значения 50 и 10, из аппроксимирующей функции

-0,9935

V = 0,2269 • S

0

вычислялось значение v, выбиралось среднее по группе значение параметра w (см. табл. 1). По заданным таким образом значениям параметров модели выполнялось моделирование развития эпидемии. В результате было получено, что только в 15% случаев относительная ошибка в предсказании максимального значения числа заболевших за неделю P не превышала 50%.

max 1

Необходимо отметить, что в работе [1] предлагается алгоритм определения значений параметров v и S0 по числу заболевших в начальный небольшой промежуток времени. Однако, как показывают результаты, представленные в таблице 2, неизбежные, даже относительно малые, ошибки в задании параметров v и S0 приводят к неприемлемым ошибкам в определении максимального числа заболевших.

По-видимому, вывод о высокой чувствительности рассмотренной математической модели к ошибкам в задании ее параметров может быть распространен и на другие, более сложные математические модели эпидемического процесса, имеющие аналогичную структуру уравнений.

Для преодоления указанных выше недостатков была выполнена следующая работа.

Заметим, что в используемой модели (1) параметр v входит только в состав произведения v  S(t). Значит, можно попытаться найти эту функцию для каждого города (для каждого значения S0) и в дальнейшем при расчетах использовать ее вместо задания отдельных значений параметров v и S0. Для этого: для каждого города, зная развитие эпидемии при значениях параметров w и v, при которых достигается наилучшее совпадение рассчитанных и наблюдавшихся значений характеристик эпидемии, с помощью второго уравнения модели (1), используя только значения числа инфицированных In = I(tn), вычисляли дискретную функцию f (tn) = v • S (tn) и ее линейную аппроксимацию f (t) = A - B • t , называемую в дальнейшем трендом. Далее все тренды, а следовательно, и города были собраны в шесть групп в зависимости от значения параметра р. В таблице 1 представлены данные только по первой группе. В шестой группе среднее значение параметра в = 1/w = 2,6, что не соответствует стандартной продолжительности за-

болевания гриппом, однако формально мы оставили эти города для дальнейшего анализа. На рисунке 3 приведены графики средних по группе трендов. На данный момент нам не удалось объяснить, почему города с различной численностью населения и расположенные в разных регионах России попали в одну группу (см. табл. 1). Их объединяет только значение параметра р. Далее было проведено моделирование развития эпидемии. Для каждого города выбрали соответствующий тренд, значения Б0, 10, V и решили второе уравнение из (1). Результаты этих расчетов представлены на рисунке 4 точками, координаты которых - ошибки расчета. По аналогии с работой [1] будем считать результаты удовлетворительными, если относительная ошибка в определении Ртх находится в интервале ± 50%, а ошибка в определении момента времени Т(Ртах), когда наступает пик эпидемии, - в интервале ± 1 неделя. Из рисунка 4 ясно, что АРтах находится в интервале ± 50% в 90% случаев, а ошибка в определении времени наступления пика эпидемии в интервале ± 1 неделя - в 78% случаев. Одновременно по обоим показателям результаты были удовлетворительными в 74% случаев (выделено черным прямоугольником). Таким образом, введение функции ^) = V  S(t) разрывает нелинейную связь в исходной системе уравнений, приводящую к их очень большой чувствительности к ошибкам в задании входных данных. Кроме того, отпадает необходимость задавать малоизвестный параметр V.

Выводы

1. Математическая модель эпидемического процесса очень чувствительна к неизбежным ошибкам в задании значений ее параметров.

2. Предложен метод математического моделирования эпидемий гриппа с использованием аппроксимирующей функции ОД = V  S(t), позволяющий преодолеть эти недостатки модели. Однако он нуждается в дальнейшей проверке на большем количестве эпидемий.

3. Для получения более точной информации о параметрах модели совершенно необходимо для каждого города сделать статистический анализ 20 - 30 эпидемий. Это позволит выполнить настройку математической модели и определить значения параметров модели для конкретного города. Таким образом удастся учесть все (или почти все) факторы, характерные для каждого города и влияющие на развитие эпидемии. Ш

Литература

- 64 с.

1. Бароян О.В., Рвачев Л.А. Математика и эпидемиология. - М.: Знание, 1977.

2. Бейли Н. Математика в биологии и медицине. - М.: Мир, 1970. - 326 с.

3. Бородулин А.И., Десятков Б.М., Шабанов А.Н. и др. Статистический разброс числа инфицируемых индивидуумов при сезонных вспышках гриппа в г. Новосибирске 2001 - 2005 годов / Проблемы охраны здоровья населения и обеспечения гигиенической и эпидемиологической безопасности окружающей среды: Сборник статей, посвященных 85-летию службы. - Новосибирск, 2007. С. 386 - 390.

4. Бородулин А.И., Десятков Б.М., Шабанов А.Н., Ярыгин А.А. Определение первых и вторых моментов в модели эпидемического процесса Барояна-Рвачева // Сибирский журнал индустриальной математики. 2007. Т. 10. № 3 (31). С. 13 - 17.

5. Бородулин А.И., Десятков Б.М., Ярыгин А.А. Статистическая модель эпидемического процесса // Сибирский журнал индустриальной математики. 2007. Т. 10. № 2 (30). С. 23 - 30.

6. Brauer F., Castillo-Chavez C. Mathematical models in population biology and epidemiology (texts in applied mathematics). - Springer, 2011. - 531 p.

7. Mode C.J., Sleeman C.K. Stochastic processes in epidemiology: HIV/AIDS, other infectious diseases and computers. - World Scientific Pub Co Inc., 2000. - 768 p.