Математическая модель эпидемии лихорадки Эбола

Яушева Ольга Александровна

Математическая модель эпидемии лихорадки Эбола

Данная работа посвящена изучению свойств математических моделей распространения инфекционных заболеваний применительно к вирусу Эбола. На основании существующих моделей предложена модель с учетом рождаемости и смертности, позволяющая учитывать возможность наличия естественного иммунитета у отдельных индивидуумов. Произведена численная реализация всех рассматриваемых моделей в системе MATLAB с учетом данных, зафиксированных при эпидемиях в Сьерра-Леоне и Либерии. По результатам моделирования выполнено сравнение предложенной модели с уже существующими.

Общественные науки в целом

Дипломы

Вуз: Санкт-Петербургский государственный университет (СПбГУ)

ID: 587d36545f1be77c40d58cd2

UUID: c643ba6e-a055-41d7-b65d-7cd8486a4289

Язык: Русский

Опубликовано: почти 8 лет назад

Просмотры: 1980

Яушева Ольга Александровна

Источник: Санкт-Петербургский государственный университет

Комментировать 0

Рецензировать 0

Скачать - 3107882 bytes

Поделиться работой

Санкт-Петербургский государственный университет Кафедра диагностики функциональных систем Яушева Ольга Александровна Выпускная квалификационная работа бакалавра Математическая модель эпидемии лихорадки Эбола Направление 010400 Прикладная математика и информатика Научный руководитель, доктор мед. наук, профессор Шишкин В.И. Научный руководитель, доктор биол. наук, профессор Кудрявцева Г.В. Санкт-Петербург 2016

Содержание Введение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . Обзор литературы . . . . . . . . . . . . . . . . . . . . . . . . . . Глава 1. Математические модели распространения заболеваний 1.1. Классическая SIR-модель . . . . . . . . . . . . . . . . . 1.2. SEIRD-модель . . . . . . . . . . . . . . . . . . . . . . . 1.3. SEIHFR-модель . . . . . . . . . . . . . . . . . . . . . . . 1.4. Построение математической модели . . . . . . . . . . . Глава 2. Численное моделирование . . . . . . . . . . . . . . . . 2.1. SIR-модель . . . . . . . . . . . . . . . . . . . . . . . . . 2.2. SEIRD-модель . . . . . . . . . . . . . . . . . . . . . . . 2.3. SEIHFR-модель . . . . . . . . . . . . . . . . . . . . . . . 2.4. SEIHFDR-модель с учетом рождаемости и смертности 2.5. Сравнение математических моделей . . . . . . . . . . . Выводы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Заключение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Список литературы . . . . . . . . . . . . . . . . . . . . . . . . . Приложение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 4 5 7 7 9 11 14 18 18 21 23 25 27 33 34 35 36

Введение Геморрагическая лихорадка Эбола является редким и опасным заболеванием, вызываемое штаммом вируса Эбола. В настоящее время вспышка эпидемии в 2014 году в Западной Африке является крупнейшей и сложнейшей с момента открытия вируса в 1976 году в Заире. Вирус вызывает острое тяжелое заболевание, которое часто заканчивается летальным исходом. По данным Всемирной Организации Здравоохранения коэффициент смертности может достигать 90% [1]. Вакцины против лихорадки Эбола на данный момент не существует, для лечения используют симптоматическую терапию. Математическое моделирование заболеваний является мощным инструментом для изучения механизмов, с помощью которых распространяется заболевание. Эпидемиологические модели служат основой для прогнозирования и оценки динамики распространения заболевания. Для сдерживания и контроля эпидемии, важно рассматривать качественные и адекватные математические модели эпидемии. В настоящее время благодаря достижениям в области математического моделирования это является реализуемой задачей. Для построения математической модели необходимо рассмотреть процесс протекания болезни. Инкубационный период, то есть период от заражения до появления первых симптомов, составляет от 2 до 21 дня. В этот период человек не заразен. Эбола распространяется при прямом контакте от человека к человеку с кровью, выделениями и другими жидкостями инфицированного человека. Прямой контакт означает, что в глаза, нос, рот, открытые раны, ссадины здорового человека попали жидкости инфицированного человека. Заболевание не передается воздушно-капельным путем. Также во время погребальных обрядов, люди имеющие прямой контакт с умершим могут быть инфицированы, так как человек остается заразным до тех пор, пока вирус находится в организме [2]. Существуют различные математические модели эпидемий. В работе, в основном будет рассмотрена SIR-модель, её модификации и дополнения. Также будут использоваться имеющиеся данные из Либерии, Сьерра-Леоне для того, чтобы параметризовать математическую модель эпидемии лихорадки Эбола. 3

Постановка задачи Целью настоящей работы являются: ∙ Описание и анализ математических моделей распространения инфекционных заболеваний. ∙ Построение математической модели распространения эпидемии лихорадки Эбола, на основе рассматриваемых моделей. ∙ Численное моделирование математических моделей распространения заболевания с реальными статистическими данными. ∙ Сравнение рассматриваемых в данной работе математических моделей. 4

Обзор литературы Теория Kermack-McKendrick представляет собой гипотезу о распространении инфекционного заболевания среди населения. В 1927 году W. O. Kermack и A. G. McKendrick опубликовали свою теорию в статье [3] и она стала источником SIR-модели. В данной модели население считается фиксированным и разделенным на три группы: здоровые или восприимчивые (от англ. susceptible) индивидуумы, способные заразиться через контакт с инфицированными (от англ. infected), а также выздоровевшие, которые перестали распространять болезнь (от англ. recovered). Данная модель до сих пор не потеряла актуальности и хорошо подходит для моделирования инфекционных заболеваний. В книге Mathematical Models in Biology [4, с. 242-254] достаточно подробно рассматривается SIR-модель, а также её модификации, такие как: SI-модель, SIS-модель, SIRS-модель. Особенностями моделей является то, что индивидуум после выздоровления может вновь попасть в класс восприимчивых индивидуумов и может быть подвержен повторному инфицированию. В статье [5] Herbert W. Hethcote были проанализированы многие математические модели распространения инфекционных заболеваний в популяциях и применены к конкретным заболеваниям. Модели рассматриваемые в данной статье основываются на классической SIR-модели. Большое внимание уделяется MSEIR модели с учетом рождаемости и смертности. В MSEIR-модели популяция делится на 5 классов: с пассивным иммунитетом, восприимчивые, латентные, инфицированные, невосприимчивые. Если мать была заражена, то некоторые антитела передаются через плаценту, таким образом новорожденный имеет временный пассивный иммунитет к инфекции. В классе 𝑀 находятся дети с пассивным иммунитетом. После того, как материнские антитела исчезают из организма новорожденного, ребенок переходит в класс восприимчивых индивидуумов. В связи с тем, что вспышка эпидемии лихорадки Эбола в Западной Африке в 2014 году считается крупнейшей, литературы на данную тему еще не так много. В 2006 году учеными J. Legrand, R.F. Grais, P.Y. Boelle, A.J. Valleron, A.Flahault была предложена [6] математическая модель распространения лихорадки Эбола среди населения основанная на SIR-модели. Особенностью данной модели является разделение населения на шесть классов: восприимчивые (от англ. susceptible), латентные (от англ. exposed), инфицированные (от англ. infected), госпитализированные (от англ. hospitalized), умершие и непохороненные (от англ. funeral), невоспри5

имчивые или иммунные (от англ. removed). Данная модель хорошо подходит для моделирования распростронения болезни, вызванной вирусом Эбола, так как достаточно полно описывает течение заболевания. Так как модель была предложена в 2006 году, до нынешней вспышки в Западной Африке, использующиеся данные для параметризации модели являются устаревшими. Эпидемиология заболевания, а также статистические данные, использующиеся в данной работе, основываются на информации опубликованной на интернет сайтах Всемирной Организации Здравоохранения [7] и Центров по контролю и профилактике заболеваний [8] в разделах посвященных заболеванию, вызванному вирусом Эбола. 6

Глава 1. Математические модели распространения заболеваний Распространение инфекционных заболеваний представляет собой сложное явление с множеством взаимодействующих факторов. Ключевая роль математической эпидемиологии заключается в создании моделей распространения патогенов. Эти модели служат в качестве математической основы для понимания сложной динамики распространения заболевания. В данной главе будут описаны и исследованы классическая SIR-модель и её модификации. 1.1. Классическая SIR-модель В классической SIR-модели (от англ. Susceptible-Infected-Removed) популяция делится на три класса: восприимчивые 𝑆(𝑡), инфицированные 𝐼(𝑡), невосприимчивые 𝑅(𝑡): ∙ 𝑆(𝑡) используется для обозначения неинфицированных индивидуумов или предрасположенных к заболеванию ∙ 𝐼(𝑡) используется для обозначения инфицированных индивидуумов, способных распространить заболевание ∙ 𝑅(𝑡) используется для обозначения индивидуумов, которые были инфицированы и выбыли из класса инфицированных в результате выздоровления или смерти Популяция считается фиксированной, таким образом 𝑆(𝑡) + 𝐼(𝑡) + 𝑅(𝑡) = 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡 = 𝑁 . SIR-модель без учета рождаемости и смертности может быть выражена следующим набором обыкновенных дифференциальных уравнений: 𝑑𝑆(𝑡) 𝛽𝑆(𝑡)𝐼(𝑡) =− 𝑑𝑡 𝑁 𝑑𝐼(𝑡) 𝛽𝑆(𝑡)𝐼(𝑡) = − 𝛾𝐼(𝑡) 𝑑𝑡 𝑁 𝑑𝑅(𝑡) = 𝛾𝐼(𝑡) 𝑑𝑡 (1) (2) (3) где 𝛽 - коэффициент, который можно интерпретировать как скорость контакта, учитывающий вероятность получения болезни в случае контакта восприимчивого индивидуума с инфицированным; 𝛾 = 1/𝑇 , где 𝑇 - время 7

болезни, коэффициент можно интерпретировать как скорость выздоровления. Начальные условия в момент времени 𝑡 = 0: 𝑆(0) = 𝑆0 > 0, 𝐼(0) = 𝐼0 > 0, 𝑅(0) = 𝑅0 > 0 Правая часть уравнения (1) описывает уменьшение популяции восприимчивых индивидуумов за счет заражения инфицированными индивидуумами восприимчивых. Первое слагаемое правой части уравнения (2) описывает увеличение популяции инфицированных индивидуумов, за счет заражения восприимчивых; второе слагаемое правой части уравнения (2) описывает уменьшение популяции инфицированных индивидуумов за счет выздоровления или смерти индивидуумов. Правая часть уравнения (3) описывает увеличение популяции невосприимчивых индивидуумом за счет выздоровления или смерти инфицированных. Наглядно схема отображена на рис. 1. Стрелки указывают возможные переходы. Параметры отражены в таблице 1. Рис. 1: Общая схема перехода индивидуума из класса в класс. инфицированные, 𝑅 - невосприимчивые. № 1 2 Переход (𝑆, 𝐼) → (𝑆 − 1, 𝐼 + 1) (𝐼, 𝑅) → (𝐼 − 1, 𝑅 + 1) 𝑆 - восприимчивые, 𝐼 - Скорость перехода (𝛽𝑆𝐼)/𝑁 𝛾𝐼 Таблица 1: Таблица переходов SIR-модели. Стоит отметить, что правые части уравнений (1), (2), (3) в сумме дают ноль, так что общий размер популяции остается неизменным. Это важное свойство данной модели. Для его сохранения, в класс невосприимчивых попадают и выздоровевшие индивидуумы и умершие. Это логичный прием, так как в обоих этих случаях индивидуум не может заразить остальную популяцию. Система (1)-(3) не является линейной и не разрешима аналитически. В настоящей работе будут рассмотрены некоторые численные методы для решения данной системы дифференциальных уравнений. Также в классическую SIR-модель представленную в параграфе мож8

но включить случай притока новорожденных в класс восприимчивых индивидуумов и оттока умерших из классов восприимчивых, инфицированных и невосприимчивых. Наглядно процесс представлен на рис.2. Рис. 2: Общая схема SIR-модели с учетом рождаемости и смертности. Положим, что 1/𝜇 - средняя продолжительность жизни. Модель можно представить набором следующих дифференциальных уравнений: 𝑑𝑆(𝑡) 𝛽𝑆(𝑡)𝐼(𝑡) =− + 𝜇𝑁 − 𝜇𝑆(𝑡) 𝑑𝑡 𝑁 𝑑𝐼(𝑡) 𝛽𝑆(𝑡)𝐼(𝑡) = − 𝛾𝐼(𝑡) − 𝜇𝐼(𝑡) 𝑑𝑡 𝑁 𝑑𝑅(𝑡) = 𝛾𝐼(𝑡) − 𝜇𝑅(𝑡) 𝑑𝑡 C начальными условиями в момент времени 𝑡 = 0: 𝑆(0) = 𝑆0 > 0, 𝐼(0) = 𝐼0 > 0, 𝑅(0) = 𝑅0 > 0 1.2. SEIRD-модель Инфекционные заболевания могут иметь скрытый или латентный период заболевания. Это период от момента, когда человек уже заражен до проявления первых симптомов. В случае лихорадки Эбола инкубационный период имеет продолжительность от 2 до 21 дня, в течение которых человек не заразен. В SEIRD-модели рассматривается ситуация, когда человек уже инфицирован, но заболевание находится в латентном периоде. Данные индивидуумы будут находится в классе латентных (от англ. exposed). В SEIRD-модели популяция делится на пять классов: восприимчивые 𝑆(𝑡), латентные 𝐸(𝑡), инфицированные 𝐼(𝑡), невосприимчивые 𝑅(𝑡), умершие 𝐷(𝑡). Классы 𝑆(𝑡), 𝐼(𝑡) определяются аналогично параграфу 1.1., 𝐸(𝑡), 𝑅(𝑡), 𝐷(𝑡) определяются следующим образом: 9

∙ 𝐸(𝑡) используется для обозначения индивидуумов, заболевание у кото- рых находится в инкубационном периоде ∙ 𝑅(𝑡) используется для обозначения выздоровевших индивидуумов ∙ 𝐷(𝑡) используется для обозначения умерших индивидуумов Популяция считается фиксированной во время вспышки, поэтому в любой момент 𝑡 популяция равна 𝑁 , т.е. 𝑆(𝑡) + 𝐸(𝑡) + 𝐼(𝑡) + 𝑅(𝑡) + 𝐷(𝑡) = 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡 = 𝑁 . SEIRD-модель может быть выражена следующим набором дифференциальных уравнений: 𝑑𝑆(𝑡) 𝛽𝑆(𝑡)𝐼(𝑡) =− 𝑑𝑡 𝑁 𝑑𝐸(𝑡) 𝛽𝑆(𝑡)𝐼(𝑡) = − 𝛿𝐸(𝑡) 𝑑𝑡 𝑁 𝑑𝐼(𝑡) = 𝛿𝐸(𝑡) − 𝛾𝐼(𝑡) − 𝜇𝐼(𝑡) 𝑑𝑡 𝑑𝑅(𝑡) = 𝛾𝐼(𝑡) 𝑑𝑡 𝑑𝐷(𝑡) = 𝜇𝐼(𝑡) 𝑑𝑡 где коэффициенты 𝛽, 𝛾 определяются аналогично параграфу 1.1.; 1/𝛿 среднее время продолжительности латентного периода; 𝜇 - коэффициент смертности. Начальные данные в момент времени 𝑡 = 0: 𝑆(0) = 𝑆0 > 0, 𝐸(0) = 𝐸0 > 0, 𝐼(0) = 𝐼0 > 0, 𝑅(0) = 𝑅0 > 0 Структура данной модели может быть представлена в виде блоксхемы рис. 3. Параметры отражены в таблице 2. Рис. 3: Общая схема SEIRD-модели: 𝑆 - восприимчивые, ванные, 𝑅 - невосприимчивые, 𝐷 - умершие. 10 𝐸 - латентные, 𝐼 - инфициро-

№ 1 2 3 4 Переход (𝑆, 𝐸) → (𝑆 − 1, 𝐸 + 1) (𝐸, 𝐼) → (𝐸 − 1, 𝐼 + 1) (𝐼, 𝑅) → (𝐼 − 1, 𝑅 + 1) (𝐼, 𝐷) → (𝐼 − 1, 𝐷 + 1) Скорость перехода (𝛽𝑆𝐼)/𝑁 𝛿𝐸 𝛾𝐼 𝜇𝐼 Таблица 2: Таблица переходов SEIRD-модели 1.3. SEIHFR-модель SEIHFR-модель делит популяцию на шесть классов: 𝑆 - восприимчивые (от англ. susceptible), 𝐸 - латентные (от англ. exposed), 𝐼 - инфицированные (от англ. infected), 𝐻 - госпитализированные (от англ. hospitalized), 𝐹 - непохороненные (от англ. funeral), 𝑅 - невосприимчивые (от англ. removed). Данная модель была предложена учеными из университета Пьера и Марии Кюри и достаточно точно описывает распространение эпидемии лихорадки Эбола, так как в данную модель включены все группы людей, участвующие в эпидемии. Остановимся подробнее на описании классов SEIHFR-модели: ∙ 𝑆(𝑡) используется для обозначения числа восприимчивых индивидуумов, находящихся в группе риска, в момент времени 𝑡 ∙ 𝐸(𝑡) используется для обозначения числа индивидуумов, заболевание которых находится в инкубационном периоде, в момент времени 𝑡 ∙ 𝐼(𝑡) используется для обозначения числа инфицированных индивидуумов, способных распространить заболевание, в момент времени 𝑡 ∙ 𝐻(𝑡) используется для обозначения числа индивидуумов, которые были госпитализированы, в момент времени 𝑡 ∙ 𝐹 (𝑡) используется для обозначения числа индивидуумов, которые умерли, но еще не похоронены, в момент времени 𝑡 ∙ 𝑅(𝑡) используется для обозначения индивидуумов, выбывших из предыдущего класса в результате выздоровления или смерти, в момент времени 𝑡 Отметим важность рассмотрения класса F в случае моделирования распространения лихорадки Эбола. Один из факторов быстрого распространения заболевания в Западной Африке это местные погребальные обряды, при которых люди имеют прямой контакт с телом умершего и могут стать переносчиком вируса Эбола. 11

Восприимчивый индивидуум 𝑆 может стать латентным 𝐸 после контакта с инфицированным 𝐼 и в дальнейшем перейти к классу инфицированных 𝐼 после инкубационного периода болезни. Часть инфицированных 𝐼 людей может быть госпитализирована 𝐻 . Для лиц из класса инфицированных 𝐼 и госпитализированных 𝐻 возможны два исхода: индивидуумы могут умереть, с шансом заразить других во время похорон 𝐹 , и перейти в класс невосприимчивых 𝑅 или они могут выздороветь и перейти в класс невосприимчивых 𝑅. SEIHFR-модель может быть выражена следующими дифференциальными уравнениями: где 𝑑𝑆(𝑡) 1 = − (𝛽𝐼 𝑆(𝑡)𝐼(𝑡) + 𝛽𝐻 𝑆(𝑡)𝐻(𝑡) + 𝛽𝐹 𝑆(𝑡)𝐹 (𝑡)) 𝑑𝑡 𝑁 𝑑𝐸(𝑡) 1 = (𝛽𝐼 𝑆(𝑡)𝐼(𝑡) + 𝛽𝐻 𝑆(𝑡)𝐻(𝑡) + 𝛽𝐹 𝑆(𝑡)𝐹 (𝑡)) − 𝛼𝐸(𝑡) 𝑑𝑡 𝑁 𝑑𝐼(𝑡) = 𝛼𝐸(𝑡) − (𝛾𝐻 𝜃1 + 𝛾𝐼 (1 − 𝜃1 )(1 − 𝛿1 ) + 𝛾𝐷 (1 − 𝜃1 )𝛿1 )𝐼(𝑡) 𝑑𝑡 𝑑𝐻(𝑡) = 𝛾𝐻 𝜃1 𝐼(𝑡) − (𝛾𝐷𝐻 𝛿2 + 𝛾𝐼𝐻 (1 − 𝛿2 ))𝐻(𝑡) 𝑑𝑡 𝑑𝐹 (𝑡) = 𝛾𝐷 (1 − 𝜃1 )𝛿1 𝐼(𝑡) + 𝛾𝐷𝐻 𝛿2 𝐻(𝑡) − 𝛾𝐹 𝐹 (𝑡) 𝑑𝑡 𝑑𝑅(𝑡) = 𝛾𝐼 (1 − 𝜃1 )(1 − 𝛿1 )𝐼(𝑡) + 𝛾𝐼𝐻 (1 − 𝛿2 )𝐻(𝑡) + 𝛾𝐹 𝐹 (𝑡) 𝑑𝑡 - коэффициент контакта в сообществе; ∙ 𝛽𝐻 - коэффициент контакта в госпитале; ∙ 𝛽𝐹 - коэффициент контакта во время похорон; ∙ 1/𝛼 - средняя продолжительность инкубационного периода; ∙ 1/𝛾𝐻 - средняя продолжительность периода от появления первых симптомов до госпитализации; ∙ 1/𝛾𝐷𝐻 - средняя продолжительность периода от госпитализации до смерти; ∙ 1/𝛾𝐼 - средняя продолжительность инфекционного периода для выздоровевшего; ∙ 1/𝛾𝐷 - средняя продолжительность инфекционного периода для умершего; ∙ 𝛽𝐼 12

- средняя продолжительность периода от госпитализации до выздоровления; ∙ 1/𝛾𝐹 - средняя продолжительность периода от смерти до погребения; ∙ 𝜃 - доля случаев госпитализации; ∙ 𝛿 - коэффициент летальности. ∙ 1/𝛾𝐼𝐻 Коэффициенты 𝜃1, 𝛿1, 𝛿2 вычисляются следующим образом: 𝜃1 = 𝛿1 = 𝜃(𝛾𝐼 (1 − 𝛿1 ) + 𝛾𝐷 𝛿1 ) 𝜃(𝛾𝐼 (1 − 𝛿1 ) + 𝛾𝐷 𝛿1 ) + (1 − 𝜃)𝛾𝐻 𝛿𝛾𝐼 𝛿𝛾𝐼𝐻 , 𝛿2 = 𝛿𝛾𝐼 + (1 − 𝛿)𝛾𝐷 𝛿𝛾𝐼𝐻 + (1 − 𝛿)𝛾𝐷𝐻 В момент времени 𝑡 = 0 начальные условия выглядят следующим образом: 𝑆(0) = 𝑆0 > 0, 𝐸(0) = 𝐸0 > 0, 𝐼(0) = 𝐼0 > 0, 𝐻(0) = 𝐻0 > 0, 𝐹 (0) = 𝐹0 > 0, 𝑅(0) = 𝑅0 > 0 Размер популяции 𝑁 фиксированный, т.е. 𝑆(𝑡) + 𝐸(𝑡) + 𝐼(𝑡) + 𝐻(𝑡) + 𝐹 (𝑡) + 𝑅(𝑡) = 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡 = 𝑁 . Структура модели и фазы заболевания представлены на рис. 4. Параметры отражены в таблице 3. Рис. 4: Общая схема SEIHFR-модели с указанием номера перехода: 𝑆 - восприимчивые, 𝐸 - латентные, 𝐼 - инфицированные, 𝐻 - госпитализированные, 𝐹 - непохороненные, 𝑅 - невосприимчивые. 13

№ 1 2 3 4 5 6 7 8 Переход Скорость перехода (𝑆, 𝐸) → (𝑆 − 1, 𝐸 + 1) (𝛽𝐼 𝑆𝐼 + 𝛽𝐻 𝑆𝐻 + 𝛽𝐹 𝑆𝐹 )/𝑁 (𝐸, 𝐼) → (𝐸 − 1, 𝐼 + 1) 𝛼𝐸 (𝐼, 𝐻) → (𝐼 − 1, 𝐻 + 1) 𝛾𝐻 𝜃1 𝐼 (𝐻, 𝐹 ) → (𝐻 − 1, 𝐹 + 1) 𝛾𝐷𝐻 𝛿2 𝐻 (𝐹, 𝑅) → (𝐹 − 1, 𝑅 + 1) 𝛾𝐹 𝐹 (𝐼, 𝑅) → (𝐼 − 1, 𝑅 + 1) 𝛾𝐼 (1 − 𝜃1 )(1 − 𝛿1 )𝐼 (𝐼, 𝐹 ) → (𝐼 − 1, 𝐹 + 1) 𝛿1 (1 − 𝜃1 )𝛾𝐷 𝐼 (𝐻, 𝑅) → (𝐻 − 1, 𝑅 + 1) 𝛾𝐼𝐻 (1 − 𝛿2 )𝐻 Таблица 3: Таблица переходов SEIHFR-модели 1.4. Построение математической модели На основе моделей, рассмотренных в параграфах 1.1, 1.2, 1.3 построим математическую модель распространения эпидемии лихорадки Эбола. Также оценим как влияют дополнительно рассмотренные аспекты течения болезни. Для этого подробнее остановимся на эпидемиологии заболевания. До конца не ясен вопрос, имеют ли люди переболевшие лихорадкой Эбола пожизненный иммунитет. Согласно данным Центров по контролю за заболеваниями [9] у людей, переживших заболевание, были обнаружены антитела к вирусу спустя 10 лет после выздоровления. Таким образом, в построении модели будем предполагать, что люди переболевшие лихорадкой Эбола, не могут заразиться вновь. Во время вспышки эпидемии лихорадки в 2010 году группа ученных протестировала кровь людей из случайно выбранных деревень в Габоне и обнаружила, что антитела к вирусу Эбола присутствовали в более чем 15% образцов и этот показатель увеличился до почти 20% в лесных районах страны. Участники исследования не имели контакт с инфицированными людьми из-за чего ученные предположили, что люди выработали антитела к вирусу Эбола после контакта с вирусом путем употребления фруктов, загрязненных слюной инфицированных летучих мышей [10]. При построении математической модели будем предполагать, что часть населения имеет естественный иммунитет к заболеванию. Разделим популяцию на семь классов: восприимчивые 𝑆 , латентные 𝐸 , инфицированные 𝐼 , госпитализированные 𝐻 , умершие и непохороненные 𝐹 , умершие и похороненные 𝐷, иммунные 𝑅. Определим классы 𝑆(𝑡), 𝐸(𝑡), 𝐼(𝑡), 𝐻(𝑡), 𝐹 (𝑡) аналогично параграфу 1.3, а классы 𝐷(𝑡), 𝑅(𝑡) следующим образом: ∙ 𝐷(𝑡) будем использовать для обозначения числа умерших и похоронен14

ных индивидуумов, которые перенесли заболевание, в момент времени 𝑡; ∙ 𝑅(𝑡) будем использовать для обозначения числа индивидуумов, которые перенесли заболевание, а также индивидуумов с естественным иммунитетом в момент времени 𝑡. Взаимодействие между классами будет происходить следующим образом: восприимчивый индивидуум 𝑆 может перейти в класс латентных 𝐸 после контакта с инфицированным 𝐼 в случае заражения. Индивидуумы из класса инфицированных 𝐼 имеют три варианта дальнейшего движения по цепочке: переход в класс непохороненных 𝐹 , затем в класс похороненных 𝐷; переход в класс госпитализированных 𝐻 , из которого в случае выздоровления переход в класс 𝑅, в случае смерти в класс 𝐹 , затем в класс 𝐷; переход в класс иммунных 𝑅. Построим модель с учетом рождаемости и смертности. Положим, 1/𝜇 - средняя продолжительность жизни. Рассмотрим коэффициенты 𝜇1 = 𝑘𝜇, 𝜇2 = (1 − 𝑘)𝜇, где 𝑘 - доля индивидуумов с естественным иммунитетом. Тогда классы восприимчивых 𝑆 и иммунных 𝑅 индивидуумов могут пополняться за счет рождаемости (отметим, что классы умерших и непохороненных 𝐹 и умерших и похороненных 𝐷 не могу на это влиять). Смертность будем рассматривать в классах восприимчивых 𝑆 , латентных 𝐸 , инфицированных 𝐼 , госпитализированных 𝐻 , иммунных 𝑅. В момент времени 𝑡 = 0 часть индивидуумов равная 𝑘𝑁 будет находиться в классе иммунных 𝑅, т.е. в момент начала вспышки некоторая часть населения имеет естественный иммунитет и не может заразиться при контакте с инфицированными. Рис. 5: Общая схема модели учетом рождаемости и смертности: 𝑆 - восприимчивые, 𝐸 - латентные, 𝐼 - инфицированные, 𝐻 - госпитализированные, 𝐹 - умершие и непохороненные , 𝐷 - умершие и похороненные, 𝑅 - иммунные. 15

Модель можно представить в виде блок-схемы, которая изображена на рис. 5. Параметры отражены в таблице 4. № Переход Скорость перехода 1 (𝑆, 𝐸) → (𝑆 − 1, 𝐸 + 1) (𝛽𝐼 𝑆𝐼 + 𝛽𝐻 𝑆𝐻 + 𝛽𝐹 𝑆𝐹 )/𝑁 2 (𝐸, 𝐼) → (𝐸 − 1, 𝐼 + 1) 𝛼𝐸 3 (𝐼, 𝐻) → (𝐼 − 1, 𝐻 + 1) 𝛾𝐻 𝜃1 𝐼 4 (𝐻, 𝐹 ) → (𝐻 − 1, 𝐹 + 1) 𝛾𝐷𝐻 𝛿2 𝐻 5 (𝐹, 𝐷) → (𝐹 − 1, 𝐷 + 1) 𝛾𝐹 𝐹 6 (𝐻, 𝑅) → (𝐻 − 1, 𝑅 + 1) 𝛾𝐼𝐻 (1 − 𝛿2 )𝐻 7 (𝐼, 𝐹 ) → (𝐼 − 1, 𝐹 + 1) 𝛾𝐷 (1 − 𝜃1 )𝛿1 𝐼 8 (𝐼, 𝑅) → (𝐼 − 1, 𝑅 + 1) 𝛾𝐼 (1 − 𝜃1 )(1 − 𝛿1 )𝐼 Таблица 4: Таблица переходов модели Построение скоростей перехода по цепочке: 1. (𝛽𝐼 𝑆(𝑡)𝐼(𝑡))/𝑁 описывает процесс изменения размера популяции со временем за счет заражения восприимчивых индивидуумов инфицированными; (𝛽𝐻 𝑆(𝑡)𝐻(𝑡))/𝑁 описывает процесс изменения размера популяции со временем за счет заражения восприимчивых индивидуумов госпитализированными; (𝛽𝐹 𝑆(𝑡)𝐹 (𝑡))/𝑁 описывает процесс изменения размера популяции со временем за счет заражения восприимчивых индивидуумов умершими во время погребальных обрядов. 2. 𝛼𝐸(𝑡) описывает процесс изменения размера популяции со временем за счет окончания инкубационного периода заболевания у индивидуума и перехода заболевания в инфекционную стадию. 3. 𝛾𝐻 𝜃1𝐼(𝑡) описывает процесс изменения размера популяции со временем за счет госпитализации инфицированных индивидуумов. 4. 𝛾𝐷𝐻 𝛿2𝐻(𝑡) описывает процесс изменения размера популяции со временем за счет смерти госпитализированных индивидуумов. 5. 𝛾𝐹 𝐹 (𝑡) описывает процесс изменения размера популяции со временем за счет погребения умерших индивидуумов. 6. 𝛾𝐼𝐻 (1 − 𝛿2)𝐻(𝑡) описывает процесс изменения размера популяции со временем за счет выздоровления госпитализированных индивидуумов. 7. 𝛾𝐷 (1 − 𝜃1)𝛿1𝐼(𝑡) описывает процесс изменения размера популяции со временем за счет смерти инфицированного индивидуума вследствие болезни. 8. 𝛾𝐼 (1 − 𝜃1)(1 − 𝛿1)𝐼(𝑡) описывает процесс изменения размера популяции со временем за счет выздоровления инфицированных индивидуумов. Рождаемость в классах 𝑆 и 𝑅 описывается как 𝜇2(𝑁 − 𝐹 (𝑡) − 𝐷(𝑡)) и 𝜇1 (𝑁 − 𝐹 (𝑡) − 𝐷(𝑡)) соответственно. Смертность, в свою очередь, в классах 16

𝐸, 𝐼 , 𝐻 , 𝑅 описывается как 𝜇𝐸(𝑡), 𝜇𝐼(𝑡), 𝜇𝐻(𝑡), 𝜇𝑅(𝑡). Как упоминалось в параграфе 1.1, важным свойством классической SIR-модели является фиксированная популяция. При построении дифференциальных уравнений, будем следить за сохранением этого свойства. Описанную модель можно представить в виде набора следующих дифференциальных уравнений: 1 𝑑𝑆(𝑡) = − (𝛽𝐼 𝑆(𝑡)𝐼(𝑡) + 𝛽𝐻 𝑆(𝑡)𝐻(𝑡) + 𝛽𝐹 𝑆(𝑡)𝐹 (𝑡)) + 𝜇2 (𝑁 − 𝐹 (𝑡)− 𝑑𝑡 𝑁 − 𝐷(𝑡)) − 𝜇𝑆(𝑡) 𝑑𝐸(𝑡) 1 = (𝛽𝐼 𝑆(𝑡)𝐼(𝑡) + 𝛽𝐻 𝑆(𝑡)𝐻(𝑡) + 𝛽𝐹 𝑆(𝑡)𝐹 (𝑡)) − 𝛼𝐸(𝑡) − 𝜇𝐸(𝑡) 𝑑𝑡 𝑁 𝑑𝐼(𝑡) = 𝛼𝐸(𝑡) − (𝛾𝐻 𝜃1 + 𝛾𝐼 (1 − 𝜃1 )(1 − 𝛿1 ) + 𝛾𝐷 (1 − 𝜃1 )𝛿1 )𝐼(𝑡) − 𝜇𝐼(𝑡) 𝑑𝑡 𝑑𝐻(𝑡) = 𝛾𝐻 𝜃1 𝐼(𝑡) − (𝛾𝐷𝐻 𝛿2 + 𝛾𝐼𝐻 (1 − 𝛿2 ))𝐻(𝑡) − 𝜇𝐻(𝑡) 𝑑𝑡 𝑑𝐹 (𝑡) = 𝛾𝐷 (1 − 𝜃1 )𝛿1 𝐼(𝑡) + 𝛾𝐷𝐻 𝛿2 𝐻(𝑡) − 𝛾𝐹 𝐹 (𝑡) 𝑑𝑡 𝑑𝐷(𝑡) = 𝛾𝐹 𝐹 (𝑡) 𝑑𝑡 𝑑𝑅(𝑡) = 𝛾𝐼 (1 − 𝜃1 )(1 − 𝛿1 )𝐼(𝑡) + 𝛾𝐼𝐻 (1 − 𝛿2 )𝐻(𝑡) + 𝜇1 (𝑁 (𝑡) − 𝐹 (𝑡)− 𝑑𝑡 − 𝐷(𝑡)) − 𝜇𝑅(𝑡) где 𝛽𝐼 , 𝛽𝐻 , 𝛽𝐹 , 𝛼, 𝛾𝐻 , 𝛾𝐷𝐻 , 𝛾𝐼 , 𝛾𝐷 , 𝛾𝐼𝐻 , 𝛾𝐹 , 𝜃1, 𝛿1, 𝛿2 определяются аналогично параграфу 1.3. Начальные условия в момент времени 𝑡 = 0: 𝑆(0) = 𝑆0 > 0, 𝐸(0) = 𝐸0 > 0, 𝐼(0) = 𝐼0 > 0, 𝐻(0) = 𝐻0 > 0, 𝐹 (0) = 𝐹0 > 0, 𝑅(0) = 𝑘𝑁 + 𝑅0 > 0, 𝐷(0) = 𝐷0 > 0 Проверим, что размер популяция остается неизменным 𝑆(𝑡) + 𝐸(𝑡) + 𝐼(𝑡) + 𝐻(𝑡) + 𝐹 (𝑡) + 𝐷(𝑡) + 𝑅(𝑡) = 𝑁 = 𝑐𝑜𝑛𝑠𝑡, т.е. 𝑑𝑆(𝑡) 𝑑𝐸(𝑡) 𝑑𝐼(𝑡) 𝑑𝐻(𝑡) 𝑑𝐹 (𝑡) 𝑑𝐷(𝑡) 𝑑𝑅(𝑡) + + + + + + =0 𝑑𝑡 𝑑𝑡 𝑑𝑡 𝑑𝑡 𝑑𝑡 𝑑𝑡 𝑑𝑡 17

Глава 2. Численное моделирование В данном разделе будут представлены результаты численного моделирования рассматриваемых в главе 1 математических моделей. Программы реализованы с помощью пакета прикладных программ MATLAB. Будет использоваться функция ode45, которая предназначена для численного интегрирования систем обыкновенных дифференциальных уравнений, используя формулы Рунге-Кутты четвертого и пятого порядка. Часть кода программ представлена в приложении. Для параметризации математических моделей распространения эпидемии лихорадки Эбола будем использовать данные по двум странам: Либерии и Сьерра-Леоне. В таблице 5 приведен начальный набор параметров. Параметры Коэффициент контакта в обществе Коэффициент контакта в госпитале Коэффициент контакта с умершими Инкубационный период Время до госпитализации Время от госпитализации до смерти Продолжительность погребальных обрядов Инфекционный период, для выздоровевших Инфекционный период, для умерших Время от госпитализации до выздоровления Доля летальных случаев 𝛿1 𝛿2 Доля случаев госпитализации 𝜃1 Средняя продолжительность жизни Население Сьерра-Леоне 0.128 0.080 0.111 10 дней 4.12 дней 6.26 дней 4.5 дней 20 дней 10.38 дней 15.88 дней 87% 0.777 0.725 41,2% 0.197 45.6 лет 6.38 млн Либерия 0.160 0.062 0.489 12 дней 3.24 дней 10.07 дней 2.01 дней 15 дней 13.31 дней 15.88 дней 57% 0.54 0.457 51,6% 0.197 60 лет 4.45 млн Таблица 5: Таблица параметров по Сьерра-Леоне и Либерии. 2.1. SIR-модель В качестве параметров модели возьмем 𝛽 = 0.128, 𝛾 = 0.0963 и рассмотрим промежуток времени в 365 дней. Общий размер популяции равен 𝑁 = 6380000. Процесс распространения заболевания начнется в случае 𝐼(0) = 𝐼0 > 0. Рассмотрим различные начальные данные 𝑆(0) = 𝑆0 , 𝐼(0) = 𝐼0 , 𝑅(0) = 𝑅0 , где количество инфицированных будем Сьерра-Леоне. 18

брать в процентном соотношении к размеру общей популяции. Рис. 6: Графики распространения заболевания для Сьерра-Леоне в зависимости от доли зараженного населения: a) 0.1% населения инфицированы, b) 1% населения инфицированы, c) 5% населения инфицированы, d) 10% населения инфицированы. Синий цвет - класс восприимчивых, зеленый - инфицированных, красный - невосприимчивых. Ось абсцисс - время (в днях), ось ординат - размер популяции (в млн. человек). Рис. 7: График распространения заболевания для Сьерра-Леоне в случае одного инфицированного индивидуума в момент времени 𝑡 = 0. На рис. 6 изображены графики распространения лихорадки Эбола 19

в Сьерра-Леоне с различными начальными данными: a) 𝑆(0) = 0.999𝑁 , 𝐼(0) = 0.001𝑁 , 𝑅(0) = 0. b) 𝑆(0) = 0.99𝑁 , 𝐼(0) = 0.01𝑁 , 𝑅(0) = 0. c) 𝑆(0) = 0.95𝑁 , 𝐼(0) = 0.05𝑁 , 𝑅(0) = 0. d) 𝑆(0) = 0.9𝑁 , 𝐼(0) = 0.1𝑁 , 𝑅(0) = 0. На рис. 7 представлен график распространения заболевания, в случае одного инфицированного индивидуума, таким образом начальные данные 𝑆(0) = 𝑁 − 1, 𝐼(0) = 1, 𝑅(0) = 0. Промежуток времени 700 дней. 𝐼𝑚𝑎𝑥 = 211740 = 0.0332𝑁 в момент времени 𝑡 = 429.2876. Параметры модели: 𝑁 = 4450000, 𝛽 = 0.16, 𝛾 = 0.0751. Рассмотрим промежуток времени в 365 дней. Либерия. Рис. 8: Графики распространения заболевания для Либерии в зависимости от доли зараженного населения: a) 0.1% населения инфицированы, b) 1% населения инфицированы, c) 5% населения инфицированы, d) 10% населения инфицированы. Синий цвет - класс восприимчивых, зеленый - инфицированных, красный - невосприимчивых. Ось абсцисс - время (в днях), ось ординат - размер популяции (в млн. человек). На рис. 8 представлены графики распространения лихорадки Эбола в Либерии с различными начальными данными: a) 𝑆(0) = 0.999𝑁 , 𝐼(0) = 0.001𝑁 , 𝑅(0) = 0. b) 𝑆(0) = 0.99𝑁 , 𝐼(0) = 0.01𝑁 , 𝑅(0) = 0. c) 𝑆(0) = 0.95𝑁 , 𝐼(0) = 0.05𝑁 , 𝑅(0) = 0. d) 𝑆(0) = 0.9𝑁 , 𝐼(0) = 0.1𝑁 , 𝑅(0) = 0. 20

На рис. 9 представлен график распространения заболевания, в случае одного инфицированного индивидуума, таким образом начальные данные 𝑆(0) = 𝑁 − 1, 𝐼(0) = 1, 𝑅(0) = 0. Промежуток времени 365 дней. 𝐼𝑚𝑎𝑥 = 780250 = 0.1753𝑁 в момент времени 𝑡 = 180.9663. Рис. 9: График распространения заболевания для Либерии в случае одного инфицированного индивидуума в момент времени 𝑡 = 0. Как видно из представленных графиков, в случае невмешательства, колоссальная доля населения оказывается инфицирована и популяция восприимчивых индивидуумов значительно сокращается. 2.2. SEIRD-модель Рассмотрим промежуток в 125 дней. Параметры модели: ∙ Сьерра-Леоне: 𝑁 = 6380000, 𝛽 = 0.128, 𝛾 = 0.05, 𝛿 = 0.1, 𝜇 = 0.0963 ∙ Либерия: 𝑁 = 4450000, 𝛽 = 0.16, 𝛾 = 0.0667, 𝛿 = 0.0833, 𝜇 = 0.0751 На рис. 10 и рис. 11 представлены графики (ось абсцисс - время в днях, ось ординат - размер популяции в млн. человек) распространения эпидемии для Сьерра-Леоне и Либерии в зависимости от доли зараженного населения: a) 95% населения восприимчивые, 4% латентные, 1% инфицированы. Начальные данные: 𝑆(0) = 0.95𝑁 , 𝐸(0) = 0.04𝑁 , 𝐼(0) = 0.01𝑁 , 𝑅(0) = 0, 𝐷(0) = 0. b) 85% населения восприимчивые, 10% латентные, 5% инфицированы. Начальные данные: 𝑆(0) = 0.85𝑁 , 𝐸(0) = 0.1𝑁 , 𝐼(0) = 0.05𝑁 , 𝑅(0) = 0, 𝐷(0) = 0. 21

Рис. 10: Графики распространения заболевания для Сьерра-Леоне в зависимости от доли зараженного населения. Рис. 11: Графики распространения заболевания для Либерии в зависимости от доли зараженного населения. . 22

c) 80% населения восприимчивые, 15% латентные, 5% инфицированы. Начальные данные: 𝑆(0) = 0.8𝑁 , 𝐸(0) = 0.15𝑁 , 𝐼(0) = 0.05𝑁 , 𝑅(0) = 0, 𝐷(0) = 0. d) 70% населения восприимчивые, 20% латентные, 10% инфицированы. Начальные данные: 𝑆(0) = 0.7𝑁 , 𝐸(0) = 0.2𝑁 , 𝐼(0) = 0.1𝑁 , 𝑅(0) = 0, 𝐷(0) = 0. Как видно из представленных графиков, заболевание распространяется менее интенсивно и пик количества инфицированных индивидуумов происходит раньше в сравнении с SIR-моделью, это связано с особенностями класса латентных индивидуумов, таких как: продолжительность инкубационного периода, отсутствие возможности распространить заболевание. Данная модель точнее описывает процесс распространения заболевания. 2.3. SEIHFR-модель Параметры модели: ∙ Сьерра-Леоне: 𝑁 = 6380000, 𝛽𝐼 = 0.128, 𝛽𝐻 = 0.08, 𝛽𝐹 = 0.111, 𝛽𝐹 = 0.111, 𝛼 = 0.1, 𝛾𝐻 = 0.2427, 𝛾𝐷𝐻 = 0.1597, 𝛾𝐹 = 0.2222, 𝛾𝐼 = 0.05, 𝛾𝐷 = 0.0963, 𝛾𝐼 𝐻 = 0.0630, 𝜃1 = 0.197, 𝛿1 = 0.777, 𝛿2 = 0.725 ∙ Либерия:: 𝑁 = 4450000, 𝛽𝐼 = 0.160, 𝛽𝐻 = 0.062, 𝛽𝐹 = 0.489, 𝛽𝐹 = 0.111, 𝛼 = 0.0833, 𝛾𝐻 = 0.3086, 𝛾𝐷𝐻 = 0.0993, 𝛾𝐹 = 0.4975, 𝛾𝐼 = 0.0667, 𝛾𝐷 = 0.0751, 𝛾𝐼 𝐻 = 0.0630, 𝜃1 = 0.197, 𝛿1 = 0.54, 𝛿2 = 0.457 На рис. 12 и 13 представлены графики (ось абсцисс - время в днях, ось ординат - размер популяции в млн. человек) распространения заболевания для Сьерра-Леоне и Либерии: a) Начальные данные: 𝑆(0) = 0.95𝑁 , 𝐸(0) = 0.04𝑁 , 𝐼(0) = 0.01𝑁 , 𝐻(0) = 0, 𝐹 (0) = 0, 𝑅(0) = 0. b) Начальные данные: 𝑆(0) = 0.85𝑁 , 𝐸(0) = 0.1𝑁 , 𝐼(0) = 0.05𝑁 , 𝐻(0) = 0, 𝐹 (0) = 0, 𝑅(0) = 0. c) Начальные данные: 𝑆(0) = 0.8𝑁 , 𝐸(0) = 0.15𝑁 , 𝐼(0) = 0.05𝑁 , 𝐻(0) = 0, 𝐹 (0) = 0, 𝑅(0) = 0. d) Начальные данные: 𝑆(0) = 0.7𝑁 , 𝐸(0) = 0.2𝑁 , 𝐼(0) = 0.1𝑁 , 𝐻(0) = 0, 𝐹 (0) = 0, 𝑅(0) = 0. На рис. 14 представлен график распространения заболевания, в случае одного инфицированного индивидуума. Как видно из представленных графиков большая часть населения переходит в класс невосприимчивых, таким образом можно сделать вывод, 23

Рис. 12: Графики распространения заболевания для Сьерра-Леоне в зависимости от доли зараженного населения. Рис. 13: Графики распространения заболевания для Либерии в зависимости от доли зараженного населения. 24

Рис. 14: График распространения заболевания в случае одного инфицированного индивидуума в момент времени 𝑡 = 0: a) для Сьерра-Леоне b) для Либерии. что даже в случае одного инфицированного индивидуума в момент времени 𝑡 = 0 в случае невмешательства, значительная часть населения перестает быть восприимчивыми. 2.4. SEIHFDR-модель с учетом рождаемости и смертности Параметры модели определяются аналогично параграфу 2.3. Параметр 𝜇 = 0.00006 для Сьерра-Леоне, 𝜇 = 0.000046 для Либерии. Так как вопрос о естественном иммунитете еще не достаточно изучен и нет доступа к реальным данным, в качестве параметра возьмем 𝑘 = 0.05. На рис. 15 и 16 представлены графики распространения заболевания для Сьерра-Леоне и Либерии в зависимости от начальных данных: a) 𝑆(0) = 0.95𝑁 − 𝑘𝑁 , 𝐸(0) = 0.04𝑁 , 𝐼(0) = 0.01𝑁 , 𝐻(0) = 0, 𝐹 (0) = 0, 𝐷(0) = 0, 𝑅(0) = 𝑘𝑁 . b) 𝑆(0) = 0.85𝑁 − 𝑘𝑁 , 𝐸(0) = 0.1𝑁 , 𝐼(0) = 0.05𝑁 , 𝐻(0) = 0, 𝐹 (0) = 0, 𝐷(0) = 0, 𝑅(0) = 𝑘𝑁 . c) 𝑆(0) = 0.8𝑁 − 𝑘𝑁 , 𝐸(0) = 0.15𝑁 , 𝐼(0) = 0.05𝑁 , 𝐻(0) = 0, 𝐹 (0) = 0, 𝐷(0) = 0, 𝑅(0) = 𝑘𝑁 . d) 𝑆(0) = 0.7𝑁 − 𝑘𝑁 , 𝐸(0) = 0.2𝑁 , 𝐼(0) = 0.1𝑁 , 𝐻(0) = 0, 𝐹 (0) = 0, 𝐷(0) = 0, 𝑅(0) = 𝑘𝑁 . Как видно из представленных графиков в данном параграфе и параграфе 2.3 в случае невмешательства, большая часть населения не просто перестает быть восприимчивыми, а переходит в класс умерших. Таким образом, отсутствие вмешательства приведет к тому, что население значи25

Рис. 15: Графики распространения заболевания для Сьерра-Леоне в зависимости от доли зараженного населения. Ось абсцисс - время (в днях), ось ординат - размер популяции (в млн. человек). Рис. 16: Графики распространения заболевания для Либерии в зависимости от доли зараженного населения. Ось абсцисс - время (в днях), ось ординат - размер популяции (в млн. человек). тельно сократится. 26

2.5. Сравнение математических моделей Принцип действия моделей различен, поэтому в качестве начальных данных выбраны разные значения. Начальные данные: 1) SEIRD и SEIHFR модели: 95% населения восприимчивые, 4% латентные, 1% инфицированы. SEIHFDR-модель: (95-K)% населения восприимчивые, 4% латентные, 1% инфицированы, K% иммунные, при K=5, 10, 15. 2) SEIRD и SEIHFR модели: 70% населения восприимчивые, 20% латентные, 10% инфицированы. SEIHFDR-модель: (70-K)% восприимчивые, 20% латентные, 10% инфицированы, K% иммунные, при K=5, 10, 15. 3) SIR, SEIRD, SEIHFR модели: 99% населения восприимчивые, 1% инфицированы. SEIHFDR-модель: (99-K)% населения восприимчивые, 1% инфицированы, K% иммунные, при K=5, 10, 15. 4) SIR, SEIRD, SEIHFR модели: 80% населения восприимчивые, 20% инфицированы. SEIHFDR-модель: (80-K)% населения восприимчивые, 20% инфицированы, K% иммунные, при K=5, 10, 15. Критерий признания окончания вспышки эпидемии Эбола согласно ВОЗ, это условие, что прошло 42 дня с момента последнего зафиксированного случая заражения [11]. Положим, что окончание вспышки эпидемии это условие существования минимально возможного 𝑡𝑒𝑛𝑑, такого что 𝐼(𝑡) < 1 при 𝑡 ∈ [𝑡𝑒𝑛𝑑 , 𝑡𝑒𝑛𝑑 + 42]. Также рассмотрим пик количества инфицированных индивидуумов 𝐼(𝑡) в момент времени 𝑡 = 𝑡𝑚𝑎𝑥 и обозначим 𝐼𝑚𝑎𝑥 . Для SEIRD и SEIHFDR моделей рассмотрим сумму классов выздоровевших и умерших 𝑃 (𝑡) = 𝐷(𝑡) + 𝑅(𝑡), которая будет обозначать класс невосприимчивых индивидуумов. Для SIR и SEIHFR моделей обозначим класс невосприимчивых индивидуумов как 𝑃 (𝑡). В таблицах 6, 7, 8, 9 представлены результаты численного моделирования рассматриваемых математических моделей с различными начальными данными. Таблица 6: Таблица результатов численного моделирования для Сьерра-Леоне. Модель 𝑡𝑚𝑎𝑥 Начальные данные 1: SEIRD 11.28 SEIHFR 72.41 SEIHFDR, K=5 73.09 SEIHFDR, K=10 70.57 SEIHFDR, K=15 67.75 Начальные данные 2: 𝐼𝑚𝑎𝑥 𝑡𝑒𝑛𝑑 1.84% 5.68% 4.95% 4.23% 3.59% 627.62 522.51 553.02 583.81 616.27 27 𝑆(𝑡𝑒𝑛𝑑 ) 𝑅(𝑡𝑒𝑛𝑑 ) 𝐷(𝑡𝑒𝑛𝑑 ) 𝑃 (𝑡𝑒𝑛𝑑 ) 79.19% 25.14% 27.08% 28.93% 30.81% 7.10% 13.77% 17.7% 21.61% 13.69% 59.13% 53.35% 47.56% 20.8% 74.85% 72.91% 71.06% 69.18%

SEIRD SEIHFR SEIHFDR, K=5 SEIHFDR, K=10 SEIHFDR, K=15 4.35 11.81 9.59 8.400 8.360 11.01% 12.94% 12.74% 12.54% 12.37% 311.77 345.25 348.79 345.85 344.73 42.23% 15.65% 16.23% 16.5% 16.62% 19.73% 15.33% 19.55% 23.78% 38.02% 68.43% 63.94% 59.58% 57.76% 84.34% 83.76% 83.49% 83.37% Таблица 7: Таблица результатов численного моделирования для Либерии. Модель 𝑡𝑚𝑎𝑥 Начальные данные 1: SEIRD 30.63 SEIHFR 69.65 SEIHFDR, K=5 71.47 SEIHFDR, K=10 71.9 SEIHFDR, K=15 72.29 Начальные данные 2: SEIRD 4.14 SEIHFR 19.11 SEIHFDR, K=5 15.39 SEIHFDR, K=10 13.09 SEIHFDR, K=15 10.44 𝐼𝑚𝑎𝑥 𝑡𝑒𝑛𝑑 𝑆(𝑡𝑒𝑛𝑑 ) 𝑅(𝑡𝑒𝑛𝑑 ) 𝐷(𝑡𝑒𝑛𝑑 ) 𝑃 (𝑡𝑒𝑛𝑑 ) 1.93% 7.92% 6.98% 6.09% 5.26% 763.15 444.31 467.63 487.98 512.98 61.58% 12.5% 14.33% 15.65% 17.08% 18.06% 39.28% 41.51% 43.68% 20.35% 46.37% 42.83% 39.22% 38.41% 87.49% 85.66% 84.34% 82.91% 10.49% 12.87% 12.33% 11.92% 11.64% 355.61 338.26 345.59 347.34 349.27 32.79% 8.37% 9.39% 9.84% 10.23% 31.59% 41.48% 44.1% 46.75% 35.61% 49.12% 46.04% 43% 67.21% 91.62% 90.60% 90.15% 89.76% Таблица 8: Таблица результатов численного моделирования для Сьерра-Леоне. Модель 𝑡𝑚𝑎𝑥 Начальные данные 3: SIR 81.33 SEIRD 0 SEIHFR 124.2 SEIHFDR, K=5 127.7 SEIHFDR, K=10 134.9 SEIHFDR, K=15 139.3 Начальные данные 4: SIR 2.49 SEIRD 0 SEIHFR 0 SEIHFDR, K=5 0 SEIHFDR, K=10 0 𝐼𝑚𝑎𝑥 𝑡𝑒𝑛𝑑 4.1% 1% 4.75% 3.97% 3.20% 2.5% 562.73 928.63 608.82 655.34 709.18 781.43 53.12% 93.65% 27.15% 29.51% 31.89% 34.47% 20.14% 20% 20% 20% 20% 281.31 363.86 386.01 394.40 397.56 32.73% - 67.26% 53.06% 16.03% 30.89% 46.93% 18.98% - 81.01% 19.94% 14.81% 65.23% 80.05% 20.65% 18.95% 60.39% 79.34% 28 𝑆(𝑡𝑒𝑛𝑑 ) 𝑅(𝑡𝑒𝑛𝑑 ) 𝐷(𝑡𝑒𝑛𝑑 ) 𝑃 (𝑡𝑒𝑛𝑑 ) 2.16% 13.4% 17.22% 20.98% 4.17% 57.07% 50.87% 44.54% 46.87% 6.34% 72.84% 70.48% 68.1% 65.52%

SEIHFDR, K=15 0 20% 399.87 21.22% 23.11% 55.66% 78.77% Таблица 9: Таблица результатов численного моделирования для Либерии. Модель 𝑡𝑚𝑎𝑥 Начальные данные 3: SIR 54.18 SEIRD 0 SEIHFR 103.3 SEIHFDR, K=5 108.1 SEIHFDR, K=10 113.7 SEIHFDR, K=15 117.7 Начальные данные 4: SIR 13.27 SEIRD 0 SEIHFR 0 SEIHFDR, K=5 0 SEIHFDR, K=10 0 SEIHFDR, K=15 0 𝐼𝑚𝑎𝑥 𝑡𝑒𝑛𝑑 𝑆(𝑡𝑒𝑛𝑑 ) 𝑅(𝑡𝑒𝑛𝑑 ) 𝐷(𝑡𝑒𝑛𝑑 ) 𝑃 (𝑡𝑒𝑛𝑑 ) 18.03% 1% 7.29% 6.32% 5.4% 4.53% 357.76 1329.5 489.21 519.33 549.2 585.83 16.83% 73% 13.26% 15.26% 16.78% 18.46% 12.69% 38.87% 40.99% 43.04% 14.3% 45.86% 42.21% 38.48% 83.16% 26.99% 86.73% 84.73% 83.21% 81.53% 28.03% 20% 20% 20% 20% 20% 283.59 416.21 363.35 372.35 380.75 387.43 12.35% 41.2% 9.91% 11.21% 11.95% 12.67% 27.64% 40.67% 43.17% 45.67% 31.15% 48.1% 44.87% 41.64% 87.64% 58.79% 90.08% 88.78% 88.04% 87.32% Интерпретация результатов с начальными данными 1 и 2 Как видно из представленных таблиц, результаты численного моделирования SEIRD-модели разительно отличаются от результатов остальных моделей: эпидемия длится дольше, её пик наступает раньше, а количество умерших при этом значительно меньше. Это связано с тем, что инфицированные индивидуумы умирают быстрее, чем распространяют болезнь. Минусом модели является не достаточно точно описанный процесс перехода индивидуумов в класс умерших D, так как не учитываются все факторы влияющие на этот процесс. По этим причинам мы опустим SEIRD-модель из дальнейшего сравнения. При количестве иммунных индивидуумов 𝐾 = 5% от общего размера популяции, видны незначительные отличия SEIHFR-модели от построенной модели, а именно: пик эпидемии, в большинстве случаев, начинается раньше, при этом с меньшим количеством инфицированных; вспышка заканчивается позже, но с большим количеством восприимчивых и меньшим количеством невосприимчивых. Это связано с тем, что класс невосприимчивых содержит в себе в том числе и умерших индивидуумов, число которых, как видно из таблиц, убывает с увеличением количества иммун29

ных. Такие результаты показывают, что добавление к рассмотрению иммунных индивидуумов приводит к ожидаемым последствиям: течение эпидемии станет равномернее, а количество умерших от болезни сократится вместе с количеством инфицированных. При увеличении количества иммунных индивидуумов 𝐾 , наблюдается усиление изменений, отмеченных в предыдущем абзаце. Из этого можно сделать вывод, что в построенной нами модели увеличение доли иммунных приводит к менее опасному течению болезни, что хорошо согласуется с реальными ожиданиями. Указанные закономерности справедливы как для результатов численного моделирования лихорадки в Либерии, так и для Сьерра-Леоне. Интерпретация результатов с начальными данными 3 и 4 При таких начальных данных вновь видно заметное отличие данных SEIRD-модели. Вне зависимости от изначального количества инфицированных индивидуумов, при нулевом значении 𝐸 (латентных) эпидемия находится в своем пике уже на старте. Связано это с тем, что при таких начальных данных, скорости изменения количества людей в классах I и E являются отрицательными. В дополнение к этому, класс латентных на старте находится в точке E=0, и в течение заболевания практически не покидает малую её окрестность, из-за чего процесс течения эпидемии становится затяжным и не динамичным. За счет простоты структуры SIR-модели, и при отсутствии зависимости от класса латентных, результаты численной реализации выглядят обособленно от других рассматриваемых моделей. Такая модель не учитывает множество важных факторов, что делает многие исследования на её основе несколько оторванными от действительности. К примеру, попытки найти для SIR-модели программное управление могут привести к результатам, невыполнимым в реальных условиях. Также из полученных при моделировании данных наглядно видно, что при отсутствии латентных индивидуумов, сравнительно большое количество инфицированных приводит к относительно быстрому и монотонному затуханию заболевания. При небольшом же начальном количестве людей в классе 𝐼 в момент времени 𝑡 = 0, лихорадка протекает долго, практически все время не уходя далеко от начального положения. Из этого можно сделать вывод, что класс латентных оказывает также немалое 30

влияние на дальнейший ход течения болезни. В дополнение к этому, стоит заметить, что в момент объявления начала эпидемии маловероятно полное отсутствие латентных индивидуумов, что делает этот класс важным для рассмотрения. По аналогии с начальными данными 1 и 2, на лицо аналогичная зависимость характера течения эпидемии от количества иммунных индивидуумов. При схожем в общих чертах поведении лихорадки, имеются заметные улучшения при наличии у людей иммунитета. Что позволяет убедиться в корректности дополнительно рассмотренных факторов течения заболевания. Рис. 17: Графики изменения количества восприимчивых индивидуумов в зависимости от выбранной модели: a) в Сьерра-Леоне b) в Либерии. Рис. 18: Графики изменения количества инфицированных индивидуумов в зависимости от выбранной модели: a) в Сьерра-Леоне b) в Либерии. На рисунках 17, 18 и 19 показаны графики изменения количества восприимчивых, инфицированных и невосприимчивых индивидуумов для рассматриваемых моделей с исходными данными 2’. Ось абсцисс - время (в днях), ось ординат - размер популяции (в млн. человек). На предложенных графиках наглядно видна обособленность результатов моделирования SEIRD модели. Также можно увидеть указанные ранее изменения в поведении SEIHFRD-модели при изменении количества иммунных индивидуумов, 31

Рис. 19: Графики изменения количества невосприимчивых индивидуумов в зависимости от выбранной модели: a) в Сьерра-Леоне b) в Либерии. и приблизительно схожий характер течения эпидемии как в Сьерра-Леоне, так и в Либерии. 32

Выводы На основе результатов численного моделирования построенной нами модели, и за счет сравнения её с данными других моделей, можно сделать следующие выводы: ∙ SIR-модель, из-за своей простоты, не учитывает большое количество важных параметров, и заметно расходится в результатах с другими, более сложными моделями. С одной стороны, эта модель легко поддается анализу и изучению, но многие исследования произведенные на её основе являются оторванными от действительности. ∙ SEIRD-модель также выделяется на фоне других моделей, так как в ней процесс перехода индивидуумов в класс умерших D рассмотрен лишь упрощенно, что отразилось на результатах моделирования. ∙ Данные SEIHFR-модели показали, что, несмотря на крайне тяжелый характер заболевания, наблюдается окончание эпидемии до полного вымирания населения, при этом число потерь среди населения является крайне высоким (более 60%). ∙ Данные полученные при анализе построенной нами SEIHFDR-модели с учетом рождаемости и смертности, показали, что при увеличении доли иммунных среди населения наблюдаются следующие улучшения: длительность эпидемии снижается, уменьшается количество как погибших от заболевания так и общее число заболевших. Такие результаты показывают корректность способа добавления к рассмотрению людей с естественным иммунитетом, так как присутствуют признаки улучшения в ходе течения вспышки эпидемии. ∙ Несмотря на заметные улучшения, даже при включении в рассмотрение возможности наличия иммунитета у индивидуумов, результаты все еще остаются неутешительными. Количество погибших индивидуумов хоть и уменьшается с ростом доли иммунных среди них, но все еще остается очень большим (от 40% до 60%). Количество потерь среди населения становится приемлемы только при высоком проценте иммунных, что на текущий момент является неправдоподобным. 33

Заключение В данной работе рассмотрены модели распространения инфекционных заболеваний: SIR, SEIRD, SEIHFR. Произведен анализ этих моделей, и на их основе предложена новая SEIHFDR-модель, учитывающая возможность наличия естественного иммунитета у отдельных индивидуумов. На основании реальных данных о лихорадке Эбола в Сьерра-Леоне и Либерии во время вспышки в 2014 году, модели были реализованы численно в среде MATLAB. Полученные при моделировании результаты позволили сравнить рассматриваемые модели между собой, и показать практическую применимость предлагаемой нами модели. Сравнение к тому же показало пригодность разработанной нами модели для дальнейших исследований на устойчивость или наличие программного управления. Также в работе показано, что, несмотря на большое количество потерь среди населения, на момент окончания эпидемии во всех рассматриваемых моделях имеется некоторое количество восприимчивых индивидуумов. То есть лихорадка заканчивается через определенный промежуток времени и без стороннего вмешательства. Хотя внешнее воздействие на течение эпидемии может существенно смягчить тяжелый характер течения эпидемии. 34

Список литературы [1] Frequently asked questions on Ebola virus disease. http://www.who.int/csr/disease/ebola/faq-ebola/en/ [2] Fact sheet №103, January 2016. http://www.who.int/mediacentre/factsheets/fs103/en/ [3] Kermack, W. O.; McKendrick, A. G. A Contribution to the Mathematical Theory of Epidemics // Proceedings of the Royal Society, 1927. Vol. 115, No. A771, P.700-721. [4] Edelstein-Keshet L. Mathematical Models in Biology. Society for Industrial and Applied Mathematics, 2005. [5] Herbert W. Hethcote. The Mathematics of Infectious Diseases // SIAM Review, 2000. Vol. 42, Iss. 4, P. 599-653. [6] J. Legrand, R.F. Grais, P.Y. Boelle, A.J. Valleron, A.Flahault. Understanding the dynamics of Ebola epidemics // Epidemiology and Infection, 2007. Vol. 135, Iss. 04, P. 610-621. [7] Ebola virus disease outbreak. http://www.who.int/csr/disease/ebola/en/ [8] Ebola Virus Disease. http://www.cdc.gov/vhf/ebola/ [9] Questions and answers on Transmission (Ebola Virus Disease). http://www.cdc.gov/vhf/ebola/transmission/qas.html [10] Becquart P, Wauquier N, Mahlakoiv T, Nkoghe D, Padilla C, Souris M, Ollomo B, Gonzalez JP, De Lamballerie X, Kazanji M, Leroy EM. High Prevalence of Both Humoral and Cellular Immunity to Zaire ebolavirus among Rural Populations in Gabon // PLoS ONE, February 2010. [11] Criteria for declaring the end of the Ebola outbreak in Guinea, Liberia or Sierra Leone. http://www.who.int/csr/disease/ebola/declaration-ebolaend/en/ 35

Приложение . Подпрограмма: SIR-модель function y=xsir(t,x) global N beta gamma xsir(1)=-beta*(1/N)*x(1)*x(2); xsir(2)=beta*(1/N)*x(1)*x(2)-gamma*x(2); xsir(3)=gamma*x(2); y=[xsir(1);xsir(2);xsir(3)]; Основная часть программы: global N beta gamma; %Параметры модели: N=6380000; beta=0.128; gamma=1/10.38; %Время начала и окончания симуляции time0=0; timef=1000; %Начальные данные x0=[S0, I0, R0] x0=[0.99*N 0.01*N 0]; %Вектор x=[S(t), I(t), R(t)] в момент времени t [t, x]=ode45(’xsir’, time0, timef, x0); SEIRD-модель. Подпрограмма: function y=xseird(t,x) global N beta gamma delta mu xsierd(1)=-beta*(1/N)*x(1)*x(3); xsierd(2)=beta*(1/N)*x(1)*x(3)-delta*x(2); xsierd(3)=delta*x(2)-gamma*x(3)-mu*x(3); xsierd(4)=gamma*x(3); xsierd(5)=mu*x(3); y=[xsierd(1);xsierd(2);xsierd(3);xsierd(4);xsierd(5)]; Основная часть программы: global N beta gamma delta mu; %Параметры модели: 36

N=6380000; beta=0.128; gamma=1/20; delta=1/10; mu=1/10.38; %Время начала и окончания симуляции time0=0; timef=1000; %Начальные данные x0=[S0, E0, I0, R0, D0] x0=[0.99*N 0 0.01*N 0 0]; %Вектор x=[S(t), E(t), I(t), R(t), D(t)] в момент времени t [t, x]=ode45(’xseird’, time0, timef, x0); SEIHFR-модель. Подпрограмма: function y=xseihfr(t,x) global N betaI betaH betaF alpha gammaH gammaDH gammaF gammaI gammaD gammaIH theta1 delta1 delta2 xseihfr(1)=-betaI*(1/N)*x(1)*x(3)-betaH*(1/N)*x(1)*x(4)-betaF*(1/N)*x(1)*x(5); xseihfr(2)=betaI*(1/N)*x(1)*x(3)+betaH*(1/N)*x(1)*x(4)+ +betaF*(1/N)*x(1)*x(5)-alpha*x(2); xseihfr(3)=alpha*x(2)-(gammaH*theta1+gammaI*(1-theta1)*(1-delta1)+ +gammaD*(1-theta1)*delta1)*x(3); xseihfr(4)=gammaH*theta1*x(3)-(gammaDH*delta2+gammaIH* *(1-delta2))*x(4); xseihfr(5)=gammaD*(1-theta1)*delta1*x(3)+gammaDH*delta2*x(4)-gammaF*x(5); xseihfr(6)=gammaI*(1-theta1)*(1-delta1)*x(3)+gammaIH*(1-delta2)* *x(4)+gammaF*x(5); y=[xseihfr(1);xseihfr(2);xseihfr(3);xseihfr(4);xseihfr(5); xseihfr(6)]; Основная часть программы: global N betaI betaH betaF alpha gammaH gammaDH gammaF gammaI gammaD gammaIH theta1 delta1 delta2 %Параметры модели: N=6380000; betaI=0.128; 37

betaH=0.080; betaF=0.111; alpha=1/10; gammaH=1/4.12; gammaDH=1/6.26; gammaF=1/4.50; gammaI=1/20; gammaD=1/10.38; gammaIH=1/15.88; theta1=0.197; delta1=0.777; delta2=0.725; %Время начала и окончания симуляции time0=0; timef=1000; %Начальные данные x0=[S0, E0, I0, H0, F0, R0] x0=[0.99*N 0 0.01*N 0 0 0]; %Вектор x=[S(t), E(t), I(t), H(t), F(t), R(t)] в момент времени t [t, x]=ode45(’xseihfr’, time0, timef, x0); SEIHFDR-модель. Подпрограмма: function y=xmymodel(t,x) global N betaI betaH betaF alpha gammaH gammaDH gammaF gammaI gammaD gammaIH theta1 delta1 delta2 mu mu1 mu2 xmymodel(1)=-betaI*(1/N)*x(1)*x(3)-betaH*(1/N)*x(1)*x(4)-betaF*(1/N)*x(1)*x(5)+mu2*(x(1)+x(2)+x(3)+x(4)+x(7))-mu*x(1); xmymodel(2)=betaI*(1/N)*x(1)*x(3)+betaH*(1/N)*x(1)*x(4)+ +betaF*(1/N)*x(1)*x(5)-alpha*x(2)-mu*x(2); xmymodel(3)=alpha*x(2)-(gammaH*theta1+gammaI*(1-theta1)* *(1-delta1)+gammaD*(1-theta1)*delta1)*x(3)-mu*x(3); xmymodel(4)=gammaH*theta1*x(3)-(gammaDH*delta2+gammaIH* *(1-delta2))*x(4)-mu*x(4); xmymodel(5)=gammaD*(1-theta1)*delta1*x(3)+gammaDH*delta2*x(4)-gammaF*x(5); xmymodel(6)=gammaF*x(5); xmymodel(7)=gammaI*(1-theta1)*(1-delta1)*x(3)+gammaIH*(1-delta2)* *x(4)-mu*x(7)+mu1*(x(1)+x(2)+x(3)+x(4)+x(7)); y=[xmymodel(1);xmymodel(2);xmymodel(3);xmymodel(4);xmymodel(5); xmymodel(6);xmymodel(7)] 38

Основная часть программы: global N betaI betaH betaF alpha gammaH gammaDH gammaF gammaI gammaD gammaIH theta1 delta1 delta2 mu k mu1 mu2 N=6380000; betaI=0.128; betaH=0.080; betaF=0.111; alpha=1/10; gammaH=1/4.12; gammaDH=1/6.26; gammaF=1/4.50; gammaI=1/20; gammaD=1/10.38; gammaIH=1/15.88; theta1=0.197; delta1=0.75; delta2=0.75; mu=1/(45.6*365); k=0.15; mu1=k*mu; mu2=(1-k)*mu; %Время начала и окончания симуляции time0=0; timef=1000; %Начальные данные x0=[S0, E0, I0, H0, F0, D0, R0] x0=[0.99*N-k*N 0 0.01*N 0 0 0 k*N]; %Вектор x=[S(t), E(t), I(t), H(t), F(t), D(t), R(t)] %в момент времени t [t, x]=ode45(’xmymodel’, time0, timef, x0); 39

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв