САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
ФАКУЛЬТЕТ ПРИКЛАДНОЙ МАТЕМАТИКИ – ПРОЦЕССОВ УПРАВЛЕНИЯ
КАФЕДРА УПРАВЛЕНИЯ МЕДИКО-БИОЛОГИЧЕСКИМИ СИСТЕМАМИ
Соснина Мария Борисовна
Выпускная квалификационная работа бакалавра
Модель эпидемической волны
Направление 010400.62
Прикладная математика и информатика
Заведующий кафедрой,
доктор физ.-мат. наук,
профессор
Александров А.Ю.
Научный руководитель,
доктор физ.-мат. наук,
профессор
Колесин И.Д.
Рецензент,
доктор физ.-мат. наук,
профессор
Колпак Е.П.
Санкт-Петербург
2016
Содержание
Введение………………………………………………………………………..…3
Постановка задачи………………………………………………………….….....4
Обзор литературы……………………………………………………………..….8
Глава 1. Случай двух городов…………………………………………………..10
1.1.
Поиск ограничений………………………………………………...10
1.2.
Задание начальных данных…………………………………….….10
Глава 2. Случай трёх городов…………………………………………………...11
1.1.
Поиск ограничений………………………………………………...11
1.2.
Задание начальных данных…………………………………….….12
Глава 3. Случай четырёх городов……………………………………………....14
1.1.
Поиск ограничений………………………………………………...14
1.2.
Задание начальных данных………………………………………..16
Выводы………………………………………………………………………...…18
Заключение……………………………………………………………………….19
Список литературы……………………………………………………………....20
Приложение……………………………………………………………………....21
2
Введение
Эпидемии, охватывающие земной шар, называются пандемиями.
Пандемии гриппа А возникают вследствие генетической изменчивости
вируса. Как показывают наблюдения, новые пандемические варианты вируса
гриппа А чаще рождаются в Юго-Восточной Азии и проникают в Россию
через контакты жителей Дальнего Востока и приграничных областей [1].
Дальнейшее распространение нового варианта осуществляется частично
через общение восприимчивых жителей с уже инфицированными, частично
через транспортные потоки. При определённых условиях возникает
эпидемическая волна, продвигающаяся в направлении Урала, а затем к
центру России [2]. Это движение моделировалось группой математиков во
главе с академиком О.В. Барояном [3]. Цель моделирования состояла в
получении возможности заблаговременного предупреждения городских
властей
и
санитарно-эпидемиологических
служб
о
надвигающейся
опасности. При этом сугубо практическая цель моделирования не
предусматривала
изучения
условий
сопряжения
отдельных
сторон
эпидемического процесса. В частности, не исследовался вопрос об условиях
возникновения бегущей волны и необходимой для этого взаимосвязи
транспортных потоков. Также не ставилась и задача управления движением
волны. Между тем, изучение этих вопросов необходимо для более полного
представления
о
пространственно-временных
особенностях
развития
эпидемических процессов и особенностях управления ими. В связи с этим,
было решено выполнить теоретическое исследование, результатом которого
стало бы получение нужных условий для возникновения волны. В данной
работе делается попытка разрешить поставленные вопросы в рамках
классических моделей типа SIR с дополнением их транспортными потоками
[4]. В соответствии с этим, требуется построить дискретную модель,
воспроизводящую движении волны и найти ограничения на интенсивность
транспортных потоков, обеспечивающих её возникновение.
3
Постановка задачи
Рассмотрим сначала случай одного города. Применяя системный
подход к построению простейшей модели эпидемии, ограничимся учетом
только трех групп индивидов: восприимчивые, зараженные, иммунные.
Пусть 𝑆, 𝐼, 𝑅 - численности этих групп. Составим уравнение постоянства
общей численности всех групп 𝑆 + 𝐼 + 𝑅 = 𝐻 = 𝑐𝑜𝑛𝑠𝑡. Выберем механизм
заражения "через встречу восприимчивого с зараженным”. Тогда
развитие эпидемии будет описываться уравнениями:
𝑑𝑆
= −𝛼𝑆𝐼,
𝑑𝑡
𝑑𝐼
= 𝛼𝑆𝐼 − 𝛽𝐼,
𝑑𝑡
𝑑𝑅
{ 𝑑𝑡 = 𝛽𝑅,
𝛽 = 1⁄𝑇 ,
𝑆(0) = 𝑆 0 > 0,
𝐼(0) = 𝐼 0 > 0,
𝑅(0) = 𝑅0 ≥ 0,
𝑇 − характерная длительность болезни.
Условием возникновения эпидемии в начальный момент времени t=0
является положительная производная
𝑑𝐼
𝑑𝑡
, вычисленная в этот момент
времени:
𝑑𝐼
(0) > 0 или 𝛼𝑆(0)𝐼(0) − 𝛽𝐼(0) > 0.
𝑑𝑡
Теперь
рассмотрим
территориального
случай
расположения
n
городов.
городов,
через
В
качестве
которые
модели
проходит
эпидемическая волна, берётся последовательность n точек на прямой, где i-ая
точка соответствует i–ому городу. Транспортные потоки как переносчики
инфекции связывают каждую точку со всеми другими, различаясь
интенсивностями переноса. Для упрощения анализа предположим, что
города идентичны друг другу по численности населения, транспортному
обслуживанию и внутригородскому перемешиванию населения, а прямые
контакты жителей одного города с другим, кроме как через транспорт,
отсутствуют. Перемещение жителей из города в город будем полагать
убывающим по интенсивности по мере удаления гостевого города от
4
исходного. Учитывая такой характер перемещений, зададим симметричную
матрицу интенсивностей перемещений P с убывающими по обе стороны от
диагоналей строчными элементами:
0
𝑝
P =( 21
⋮
𝑝𝑛1
𝑝12
0
⋮
𝑝𝑛2
𝑝13
𝑝23
⋮
𝑝𝑛3
⋯ 𝑝1𝑛
… 𝑝2𝑛
),
⋱
⋮
⋯ 0
𝑝12 > 𝑝13 > ⋯ > 𝑝1𝑛 > 0,
𝑝21 > 𝑝23 > ⋯ > 𝑝2𝑛 > 0,
𝑝𝑖𝑗 = 𝑝𝑗𝑖 , 𝑖 ≠ 𝑗
…
Кроме того, потребуем 𝑝𝑖𝑗 = 𝑝𝑖+1 , 𝑖, 𝑗 = 1,2 …, что эквивалентно сдвигу
(i+1)-ой строки относительно i-ой на один элемент вправо. Этим задаётся
идентичность транспортного обслуживания n городов. Пусть 𝑆𝑖 , 𝐼𝑖 , 𝑅𝑖 соответственно число восприимчивых, инфицированных и иммунных лиц в iом городе, 𝐻𝑖 - общее их число. Без учёта транспортных связей развитие
эпидемии в i-ом городе носило бы локальный характер, допуская
классическое описание через контакты восприимчивых с инфицированными:
𝑑𝑆𝑖
𝑑𝐼𝑖
𝑑𝑅𝑖
= −𝑎𝐼𝑖 𝑆𝑖 ,
= 𝑎𝐼𝑖 𝑆𝑖 − 𝛽𝐼𝑖 ,
= 𝛽𝐼𝑖 ,
𝑖 = ̅̅̅̅̅
1, 𝑛,
𝑑𝑡
𝑑𝑡
𝑑𝑡
где 𝛽 = 1⁄𝑇, T - характерная длительность болезни. Для включения в эту
модель
транспортных
предположения
потоков
относительно
необходимо
доли
принять
дополнительные
инфицированных
лиц
среди
перемещающейся части населения. Это особенно важно в случае удалённых
друг от друга городов, когда транспортные потоки выполняют основную
роль в распространении инфекции. Предположим, что структура контингента
лиц, перемещающихся из i-ого города в j-ый (𝑆𝑖𝑗 , 𝐼𝑖𝑗 , 𝑅𝑖𝑗 ), повторяет в
определённой пропорции структуру населения i-ого города:
𝑆𝑖𝑗 : 𝐼𝑖𝑗 : 𝑅𝑖𝑗 = 𝑆𝑖 : 𝐼𝑖 : 𝑅𝑖 .
Введём коэффициент пропорциональности 𝛿𝑖𝑗 (0<𝛿𝑖𝑗 <1), отражающий
отношение вместимости транспортных средств, связывающий i-ый город с j5
м, к числу жителей i-ого города. Тогда интенсивность транспортного потока
из i-ого города в j-ый с учётом 𝑆𝑖 + 𝐼𝑖 + 𝑅𝑖 = 𝐻𝑖 получит выражение:
𝑝𝑖𝑗 = 𝑆𝑖𝑗 + 𝐼𝑖𝑗 + 𝑅𝑖𝑗 = 𝛿𝑖𝑗 𝑆𝑖𝑗 + 𝛿𝑖𝑗 𝐼𝑖𝑗 + 𝛿𝑖𝑗 𝑅𝑖𝑗 = 𝛿𝑖𝑗 𝐻𝑖𝑗 .
Аналогично для
𝑝𝑗𝑖 . В этом случае количественное изменение
структуры населения i-ого и j-ого городов за счёт притока и оттока жителей
за единицу времени ∆𝑡 = 1 составит:
∆𝑆𝑖 = −𝛿𝑖𝑗 𝑆𝑖 + 𝛿𝑗𝑖 𝑆𝑗 ,
∆𝐼𝑖 = −𝛿𝑖𝑗 𝐼𝑖 + 𝛿𝑗𝑖 𝐼𝑗 ,
∆𝑆𝑗 = 𝛿𝑖𝑗 𝑆𝑖 − 𝛿𝑗𝑖 𝑆𝑗 ,
∆𝑅𝑖 = −𝛿𝑖𝑗 𝑅𝑖 + 𝛿𝑗𝑖 𝑅𝑗 ,
∆𝐼𝑗 = 𝛿𝑖𝑗 𝐼𝑖 − 𝛿𝑗𝑖 𝐼𝑗 ,
𝑆𝑖 + 𝐼𝑖 + 𝑅𝑖 = −𝛿𝑖𝑗 𝐻𝑖 + 𝛿𝑗𝑖 𝐻𝑗 ,
∆𝑅𝑖 = 𝛿𝑖𝑗 𝑅𝑖 − 𝛿𝑗𝑖 𝑅𝑗 ,
𝑆𝑗 + 𝐼𝑗 + 𝑅𝑗 = 𝛿𝑖𝑗 𝐻𝑖 − 𝛿𝑗𝑖 𝐻𝑗 .
Вместе с тем, количественные изменения эпидемиологического
характера за тот же отрезок времени ∆𝑡 = 1 составят:
∆𝑆𝑖 = −𝑎𝑖 (𝐼𝑖 + ∑ 𝛿𝑖𝑗 𝐼𝑖 ) 𝑆𝑖 , ∆𝐼𝑖 = 𝑎𝑖 (𝐼𝑖 + ∑ 𝛿𝑖𝑗 𝐼𝑖 ) 𝑆𝑖 − 𝛽𝐼𝑖 , ∆𝑅𝑖 = 𝛽𝐼𝑖 ,
𝑗≠𝑖
𝑗≠𝑖
∆𝑆𝑗 = −𝑎𝑖 (𝐼𝑗 + ∑ 𝛿𝑗𝑖 𝐼𝑗 ) 𝑆𝑗 , ∆𝐼𝑗 = 𝑎𝑖 (𝐼𝑗 + ∑ 𝛿𝑗𝑖 𝐼𝑗 ) 𝑆𝑗 − 𝛽𝐼𝑗 , ∆𝑅𝑗 = 𝛽𝐼𝑗.
𝑗≠𝑖
Объединяя
𝑗≠𝑖
транспортные
изменения
с
эпидемиологическими
и
переходя к дифференциальным уравнениям, получим следующую модель
развития пространственно-временного процесса в системе из n городов:
𝑑𝑆𝑖
𝑑𝑡
𝑑𝐼𝑖
𝑑𝑡
= −𝑎𝑖 (𝐼𝑖 + ∑𝑗≠𝑖 𝛿𝑗𝑖 𝐼𝑗 )𝑆𝑖 − ∑𝑗≠𝑖 𝛿𝑖𝑗 𝑆𝑖 + ∑𝑗≠𝑖 𝛿𝑗𝑖 𝑆𝑗 ,
= 𝑎(𝐼𝑖 + ∑𝑗≠𝑖 𝛿𝑗𝑖 𝐼𝑗 )𝑆𝑖 − ∑𝑗≠𝑖 𝛿𝑖𝑗 𝐼𝑖 + ∑𝑗≠𝑖 𝛿𝑗𝑖 𝐼𝑗 − 𝛽𝐼𝑖 ,
{
𝑑𝑅𝑖
𝑑𝑡
(1)
= 𝛽𝐼𝑖 − ∑𝑗≠𝑖 𝛿𝑖𝑗 𝑅𝑖 + ∑𝑗≠𝑖 𝛿𝑗𝑖 𝑅𝑗 .
Предположим, что в момент времени t=0 в каждом из n городов уже
имеется некоторое количество инфицированных и задано начальное
распределение инфицированных лиц по n городам:
𝑡 = 0: 𝐼𝑖 (0) > 0, 𝑖 = 1, 𝑛.
Далее предположим, что в первом городе число инфицированных лиц
больше, чем в других:
6
𝐼1 (0) > 𝐼𝑖 (0), 𝑖 = 2, 𝑛.
(2)
Будем полагать, что в первом городе в начальный момент времени
выполняется условие начала эпидемии:
𝑡 = 0:
𝑑𝐼1
𝑑𝑡
(0) > 0,
𝑑𝐼𝑖
𝑑𝑡
(0) < 0, 𝑖 = 2, 𝑛.
(3)
Требуется найти такие условия на интенсивность транспортных
потоков 𝛿𝑖𝑗 , при которых для системы (1) выполняются условия (2) и (3).
Задача состоит в нахождении ограничений для 𝛿𝑖𝑗 , обеспечивающих
выполнение условий (2) и (3): 0 < 𝛿𝑖𝑗 < 𝛿̂
𝑖𝑗 .
Найдём последовательно решение задачи (1)-(3) для случая двух
городов, а затем для случаев трёх и четырёх.
7
Обзор литературы
История количественных наблюдений за развитием эпидемий идёт с
конца 19 века. Исторически первая математическая модель эпидемии была
создана в 1927 году У. Кермаком и А. Маккендриком и описана в статье «A
Contribution to the Mathematical Theory of Epidemics». Учёные предложили
модель для одного города с населением, неизменным во времени по
численности. Ими было предложено все население рассматриваемого города
разделить на три группы: индивиды, которые восприимчивы к данной
болезни, но здоровы («восприимчивые», susceptible) — S(t); зараженные
индивиды («инфицированные», infected) — I(t); и здоровые индивиды,
обладающие иммунитетом к данной болезни («выздоровевшие», recovered)
— R(t). Это детерминированная модель, в которой механизм заражения
реализуется через встречи восприимчивых с зараженными. Множество
моделей подобного типа, но с разным числом взаимодействующих групп
исследовалось в 80-е годы.
Движение эпидемической волны по территории СССР с захватом
большого числа городов моделировалось группой математиков во главе с
академиком О.В. Барояном в работе «Моделирование и прогнозирование
эпидемий гриппа для территории СССР». Учёными была разработана
новая методология математического моделирования эпидемий. Данная
методология основана на методе научной аналогии в отображении
эпидемического процесса (процесс «переноса» возбудителя инфекции от
больных к здоровым) с процессом «переноса» материи (энергии, импульса и
др.) в уравнениях математической физики. Система уравнений, которая
описывает развитие эпидемического процесса, представляет собой систему
нелинейных уравнений в частных производных с соответствующими
начальными и граничными условиями, весьма «схожими» с уравнениями
гидродинамики.
8
С применением этой методологии в 60-70-егоды были разработаны
уникальные модели эпидемий гриппа для территории СССР, которые
составлены
на
инфекционного
основе
основных
процесса
типа
проходящих
SIR,
где
S
стадий
состояния
-восприимчивые,
I
-
инфицированные больные, R – переболевшие.
Цель
моделирования
заблаговременного
состояла
предупреждения
в
получении
городских
властей
возможности
и
санитарно-
эпидемиологических служб о надвигающейся опасности. При этом сугубо
практическая цель моделирования не предусматривала изучения условий
сопряжения отдельных сторон эпидемического процесса. В частности, не
исследовался вопрос об условиях возникновения бегущей волны и
необходимой для этого взаимосвязи транспортных потоков. Также не
ставилась и задача управления движением волны.
В настоящее время создано множество различных моделей эпидемий.
Многие из них описаны в учебном пособии Колесина И.Д. и Житковой Е.М.
«Математические модели эпидемий» [5].
9
Глава 1. Случай двух городов
1.1. Поиск ограничений
Уравнения, описывающие изменение количества инфицированных лиц
за единицу времени ∆𝑡 = 1, для случая двух городов имеют вид:
𝑑𝐼1
= 𝑎(𝐼1 + 𝛿21 𝐼2 )𝑆1 − 𝛿12 𝐼1 + 𝛿21 𝐼2 − 𝛽𝐼1 ,
𝑑𝑡
𝑑𝐼2
= 𝑎(𝐼2 + 𝛿12 𝐼1 )𝑆2 − 𝛿21 𝐼2 + 𝛿12 𝐼1 − 𝛽𝐼2 .
𝑑𝑡
Пусть для каждого из двух городов выполняется условие (2). Для
выполнения условия (3) потребуем
𝑎(𝐼1 + 𝛿21 𝐼2 )𝑆1 − 𝛿12 𝐼1 + 𝛿21 𝐼2 − 𝛽𝐼1 > 0,
𝑎(𝐼2 + 𝛿12 𝐼1 )𝑆2 − 𝛿21 𝐼2 + 𝛿12 𝐼2 − 𝛽𝐼1 < 0.
Учитывая симметричность матрицы интенсивности перемещений,
примем 𝛿12 = 𝛿21 . Из первого неравенства следует:
𝛿12 >
𝐼1 (𝛽 − 𝛼𝑆1 )
.
𝛼𝑆1 𝐼1 − 𝐼1 + 𝐼2
Из второго неравенства следует:
𝛿12 <
𝐼2 (𝛽 − 𝛼𝑆2 )
.
𝛼𝑆2 𝐼1 − 𝐼2 + 𝐼1
Объединяя, получим следующее ограничение на выбор 𝛿12 и 𝛿21 :
𝐼1 (𝛽 − 𝛼𝑆1 )
𝐼2 (𝛽 − 𝛼𝑆2 )
< 𝛿12 <
.
𝛼𝑆1 𝐼1 − 𝐼1 + 𝐼2
𝛼𝑆2 𝐼1 − 𝐼2 + 𝐼1
1.2. Задание начальных данных
Для нахождения конкретных значений для 𝛿12 и 𝛿21 зададим
следующие начальные данные, удовлетворяющие условию (2): 𝛽 = 0,1,
𝛼 = 0,00025, 𝑆1 = 900, 𝐼1 = 100, 𝑅1 = 0, 𝑆2 = 950, 𝐼2 = 50, 𝑅2 = 0.
В соответствии с выбранными данными, получим следующий интервал
для выбора 𝛿12 и 𝛿21 :
0.016 < 𝛿12 < 0.434.
10
Выберем
ограничению.
𝛿12 = 𝛿21 = 0,03,
удовлетворяющее
полученному
Модель эпидемии для случая двух городов с учётом полученных
данных была реализована в программе Matlab (см. Приложение). Результат
данного практического эксперимента наглядно представлен на Рис.1.
Рисунок 1.
На рисунке показано, что в первом городе возникает эпидемия,
достигает своего пика при t=15, а затем наступает спад. Эпидемия во втором
городе возникает практически одновременно с первым, задержка составляет
несколько дней. При этом общее количество инфицированных во втором
городе меньше, чем в первом. В первом городе наступает процесс
выздоровления, и количество инфицированных лиц уменьшается. Затем
процесс выздоровления происходит и во втором городе.
Глава 2. Случай трёх городов
2.1. Поиск ограничений
Уравнения, описывающие изменение количества инфицированных лиц
за единицу времени ∆𝑡 = 1, для случая трёх городов имеют вид:
11
𝑑𝐼1
= 𝑎(𝐼1 + 𝛿21 𝐼2 + 𝛿31 𝐼3 )𝑆1 − 𝛿12 𝐼1 − 𝛿13 𝐼1 + 𝛿21 𝐼2 + 𝛿31 𝐼3 − 𝛽𝐼1 ,
𝑑𝑡
𝑑𝐼2
= 𝑎(𝐼2 + 𝛿12 𝐼1 + 𝛿32 𝐼3 )𝑆2 − 𝛿21 𝐼2 − 𝛿23 𝐼2 + 𝛿12 𝐼1 +𝛿32 𝐼3 − 𝛽𝐼2 ,
𝑑𝑡
𝑑𝐼3
= 𝑎(𝐼3 + 𝛿13 𝐼1 + 𝛿23 𝐼2 )𝑆3 − 𝛿31 𝐼3 − 𝛿32 𝐼3 + 𝛿13 𝐼1 + 𝛿23 𝐼2 − 𝛽𝐼3 .
𝑑𝑡
Пусть выполняется условие (2). Для выполнения условия (3) потребуем
𝑎(𝐼1 + 𝛿21 𝐼2 + 𝛿31 𝐼3 )𝑆1 − 𝛿12 𝐼1 − 𝛿13 𝐼1 + 𝛿21 𝐼2 + 𝛿31 𝐼3 − 𝛽𝐼1 > 0,
𝑎(𝐼2 + 𝛿12 𝐼1 + 𝛿32 𝐼3 )𝑆2 − 𝛿21 𝐼2 − 𝛿23 𝐼2 + 𝛿12 𝐼1 +𝛿32 𝐼3 − 𝛽𝐼2 < 0,
𝑎(𝐼3 + 𝛿13 𝐼1 + 𝛿23 𝐼2 )𝑆3 − 𝛿31 𝐼3 − 𝛿32 𝐼3 + 𝛿13 𝐼1 + 𝛿23 𝐼2 − 𝛽𝐼3 < 0.
Учитывая симметричность матрицы интенсивности перемещений,
примем 𝛿13 = 𝛿31 , 𝛿23 = 𝛿32 . Из первого неравенства следует:
𝛿13 >
−𝛼𝐼1 𝑆1 − 𝛼𝐼2 𝑆1 𝛿12 + 𝛽𝐼1
.
𝛼𝑆1 𝐼3 − 𝐼1 + 𝐼3
Из второго неравенства следует:
𝛿23 <
−𝛼𝐼2 𝑆2 − 𝛼𝐼2 𝑆2 𝛿12 + 𝐼2 𝛿12 + 𝐼1 𝛿12 + 𝛽𝐼2
.
𝛼𝑆2 𝐼3 − 𝐼2 + 𝐼3
Из третьего неравенства следует:
𝛿13 <
−𝛼𝐼3 𝑆3 − 𝛼𝐼2 𝑆3 𝛿23 + 𝐼3 𝛿32 − 𝐼2 𝛿23 + 𝛽𝐼3
𝛼𝑆3 𝐼3 − 𝐼3 + 𝐼1
Объединяя, получим следующие ограничения на выбор 𝛿13 , 𝛿31 ,
𝛿23 и 𝛿32 :
−𝛼𝐼1 𝑆1 − 𝛼𝐼2 𝑆1 𝛿12 + 𝛽𝐼1
−𝛼𝐼3 𝑆3 − 𝛼𝐼2 𝑆3 𝛿23 + 𝐼3 𝛿32 − 𝐼2 𝛿23 + 𝛽𝐼3
< 𝛿13 <
,
𝛼𝑆1 𝐼3 − 𝐼1 + 𝐼3
𝛼𝑆3 𝐼3 − 𝐼3 + 𝐼1
0 < 𝛿23 <
−𝛼𝐼2 𝑆2 − 𝛼𝐼2 𝑆2 𝛿12 + 𝐼2 𝛿12 + 𝐼1 𝛿12 + 𝛽𝐼2
.
𝛼𝑆2 𝐼3 − 𝐼2 + 𝐼3
2.2. Задание начальных данных
Для нахождения конкретных значений для 𝛿13 , 𝛿31 , 𝛿23 и 𝛿32 зададим
следующие начальные данные для третьего города, удовлетворяющие
условию (2): 𝑆3 = 985, 𝐼3 = 15, 𝑅3 = 0. Используем начальные данные для
12
первых двух городов, введённые ранее, и выбранное значение для 𝛿12 =
𝛿21 = 0,03.
В соответствии с введёнными и ранее полученными данными, получим
интервал для выбора 𝛿23 и 𝛿32 :
0 < 𝛿23 < 0,891.
Принимая во внимание, что матрица перемещений – симметричная
матрица с убывающими по обе стороны от диагоналей элементами, выберем
𝛿23 = 𝛿32 = 0,01, удовлетворяющие полученному ограничению.
Аналогично получим следующий интервал для выбора 𝛿13 и 𝛿31 ,
используя при его расчёте выбранное значение для 𝛿23 = 𝛿32 = 0,01:
−0.037 < 𝛿13 < 0,628.
Учитывая ранее обозначенные ограничения на 𝛿𝑖𝑗 (0<𝛿𝑖𝑗 <1) получим:
0 < 𝛿13 < 0,628.
Так как матрица перемещений – симметричная матрица с убывающими
по обе стороны от диагоналей элементами, выберем 𝛿13 = 𝛿31 =
0,1, удовлетворяющие полученному выше ограничению.
Результат реализации модели для случая трёх городов в программе
Matlab представлен на Рис.2. На этом рисунке показано, что вспышка
эпидемии в третьем городе возникает не сразу после вспышки во втором,
требуется некоторое время для того, чтобы количество инфицированных лиц
в третьем городе возрасло. Об этой особенности говорит более плавное
возрастание графика для третьего города. Уровень вспышки в третьем городе
заметно ниже, чем в двух предыдущих. Это можно объяснить меньшим
количеством инфицированных в этом городе.
13
Рисунок 2.
Глава 3. Случай четырёх городов
3.1. Поиск ограничений
Уравнения, описывающие изменение количества инфицированных лиц
за единицу времени ∆𝑡 = 1, для случая четырёх городов имеют вид:
𝑑𝐼1
= 𝑎(𝐼1 + 𝛿21 𝐼2 + 𝛿31 𝐼3 + 𝛿41 𝐼4 )𝑆1 − 𝛿12 𝐼1 − 𝛿13 𝐼1 − 𝛿14 𝐼1 + 𝛿21 𝐼2 + 𝛿31 𝐼3
𝑑𝑡
+ 𝛿41 𝐼4 − 𝛽𝐼1 ,
𝑑𝐼2
= 𝑎(𝐼2 + 𝛿12 𝐼1 + 𝛿32 𝐼3 + 𝛿42 𝐼4 )𝑆2 − 𝛿21 𝐼2 − 𝛿23 𝐼2 − 𝛿24 𝐼2 + 𝛿12 𝐼1 + 𝛿32 𝐼3
𝑑𝑡
+ 𝛿42 𝐼4 − 𝛽𝐼2 ,
𝑑𝐼3
= 𝑎(𝐼3 + 𝛿13 𝐼1 + 𝛿23 𝐼2 + 𝛿43 𝐼4 )𝑆3 − 𝛿31 𝐼3 − 𝛿32 𝐼3 − 𝛿34 𝐼3 + 𝛿13 𝐼1 + 𝛿23 𝐼2
𝑑𝑡
+ 𝛿43 𝐼4 − 𝛽𝐼3 ,
𝑑𝐼4
= 𝑎(𝐼4 + 𝛿14 𝐼1 + 𝛿24 𝐼2 + 𝛿34 𝐼3 )𝑆4 − 𝛿41 𝐼4 − 𝛿42 𝐼4 + 𝛿43 𝐼4 + 𝛿14 𝐼1 + 𝛿24 𝐼2
𝑑𝑡
+ 𝛿34 𝐼3 − 𝛽𝐼4 .
14
Пусть выполняется условие (2). Для выполнения условия (3) потребуем
𝑎(𝐼1 + 𝛿21 𝐼2 + 𝛿31 𝐼3 + 𝛿41 𝐼4 )𝑆1 − 𝛿12 𝐼1 − 𝛿13 𝐼1 − 𝛿14 𝐼1 + 𝛿21 𝐼2 + 𝛿31 𝐼3 + 𝛿41 𝐼4
− 𝛽𝐼1 > 0,
𝑎(𝐼2 + 𝛿12 𝐼1 + 𝛿32 𝐼3 + 𝛿42 𝐼4 )𝑆2 − 𝛿21 𝐼2 − 𝛿23 𝐼2 − 𝛿24 𝐼2 + 𝛿12 𝐼1 + 𝛿32 𝐼3 +
𝛿42 𝐼4 − 𝛽𝐼2 < 0,
𝑎(𝐼3 + 𝛿13 𝐼1 + 𝛿23 𝐼2 + 𝛿43 𝐼4 )𝑆3 − 𝛿31 𝐼3 − 𝛿32 𝐼3 − 𝛿34 𝐼3 + 𝛿13 𝐼1 + 𝛿23 𝐼2 + 𝛿43 𝐼4
− 𝛽𝐼3 < 0,
𝑎(𝐼4 + 𝛿14 𝐼1 + 𝛿24 𝐼2 + 𝛿34 𝐼3 )𝑆4 − 𝛿41 𝐼4 − 𝛿42 𝐼4 + 𝛿43 𝐼4 + 𝛿14 𝐼1 + 𝛿24 𝐼2 + 𝛿34 𝐼3
− 𝛽𝐼4 < 0.
Учитывая симметричность матрицы интенсивности перемещений,
примем 𝛿14 = 𝛿41 , 𝛿24 = 𝛿42 , 𝛿34 = 𝛿43 . Из первого неравенства следует:
𝛿14 >
−𝛼𝐼1 𝑆1 − 𝛼𝐼2 𝑆1 𝛿21 − 𝛼𝐼3 𝑆1 𝛿31 + 𝐼1 𝛿12 − 𝐼1 𝛿13 + 𝐼2 𝛿21 + 𝐼3 𝛿31 + 𝛽𝐼1
.
𝛼𝑆1 𝐼4 + 𝐼4 − 𝐼1
Из второго неравенства следует:
𝛿24 <
−𝛼𝐼1 𝑆1 − 𝛼𝐼2 𝑆1 𝛿21 − 𝛼𝐼3 𝑆1 𝛿31 + 𝐼1 𝛿12 − 𝐼1 𝛿13 + 𝐼2 𝛿21 + 𝐼3 𝛿31 + 𝛽𝐼1
.
𝛼𝑆1 𝐼4 + 𝐼4 − 𝐼1
Из третьего неравенства следует:
𝛿34 <
−𝛼𝐼1 𝑆1 − 𝛼𝐼2 𝑆1 𝛿21 − 𝛼𝐼3 𝑆1 𝛿31 + 𝐼1 𝛿12 − 𝐼1 𝛿13 + 𝐼2 𝛿21 + 𝐼3 𝛿31 + 𝛽𝐼1
.
𝛼𝑆1 𝐼4 + 𝐼4 − 𝐼1
Из третьего неравенства следует:
𝛿34 <
−𝛼𝐼1 𝑆1 − 𝛼𝐼2 𝑆1 𝛿21 − 𝛼𝐼3 𝑆1 𝛿31 + 𝐼1 𝛿12 − 𝐼1 𝛿13 + 𝐼2 𝛿21 + 𝐼3 𝛿31 + 𝛽𝐼1
.
𝛼𝑆1 𝐼4 + 𝐼4 − 𝐼1
Объединяя, получим следующее ограничение на выбор 𝛿14 , 𝛿41 , 𝛿24 ,
𝛿42 , 𝛿34 и 𝛿43 :
−𝛼𝐼1 𝑆1 − 𝛼𝐼2 𝑆1 𝛿12 + 𝛽𝐼1
−𝛼𝐼3 𝑆3 − 𝛼𝐼2 𝑆3 𝛿23 + 𝐼3 𝛿32 − 𝐼2 𝛿23 + 𝛽𝐼3
< 𝛿14 <
,
𝛼𝑆1 𝐼3 − 𝐼1 + 𝐼3
𝛼𝑆3 𝐼3 − 𝐼3 + 𝐼1
𝛿24 <
−𝛼𝐼1 𝑆1 − 𝛼𝐼2 𝑆1 𝛿21 − 𝛼𝐼3 𝑆1 𝛿31 + 𝐼1 𝛿12 − 𝐼1 𝛿13 + 𝐼2 𝛿21 + 𝐼3 𝛿31 + 𝛽𝐼1
.
𝛼𝑆1 𝐼4 + 𝐼4 − 𝐼1
15
𝛿34 <
−𝛼𝐼1 𝑆1 − 𝛼𝐼2 𝑆1 𝛿21 − 𝛼𝐼3 𝑆1 𝛿31 + 𝐼1 𝛿12 − 𝐼1 𝛿13 + 𝐼2 𝛿21 + 𝐼3 𝛿31 + 𝛽𝐼1
.
𝛼𝑆1 𝐼4 + 𝐼4 − 𝐼1
3.2. Задание начальных данных
Для
нахождения
конкретных
значений
для
𝛿14 , 𝛿41 ,
𝛿24 , 𝛿42 , 𝛿43 и 𝛿34 зададим следующие начальные данные для четвёртого
города,
удовлетворяющие
условию
(2):
𝑆4 = 999, 𝐼4 = 1, 𝑅4 = 0.
Используем начальные данные для первых трёх городов, введённые ранее, и
выбранное значение для 𝛿12 = 𝛿21 = 0,03, 𝛿13 = 𝛿31 = 0,02 и 𝛿23 = 𝛿32 =
0,01.
В соответствии с введёнными и ранее полученными данными, получим
интервал для выбора 𝛿14 и 𝛿41 :
0,0014 < 𝛿14 < 0,245.
Принимая во внимание, что матрица перемещений – симметричная
матрица с убывающими по обе стороны от диагоналей элементами, выберем
𝛿14 = 𝛿41 = 0,001, удовлетворяющие полученному ограничению.
Аналогично получим следующий интервал для выбора 𝛿42 и 𝛿24 ,
используя при его расчёте выбранное значение для 𝛿14 = 𝛿41 = 0,01:
−0.019 < 𝛿42 < 0,508.
Учитывая ранее обозначенные ограничения на 𝛿𝑖𝑗 (0<𝛿𝑖𝑗 <1) получим:
0 < 𝛿42 < 0,508.
Так как матрица перемещений – симметричная матрица с убывающими
по обе стороны от диагоналей элементами, выберем 𝛿42 = 𝛿24 =
0,01, удовлетворяющие полученному выше ограничению.
Аналогично получим следующий интервал для выбора 𝛿43 и 𝛿34 :
0 < 𝛿43 < 0,812.
Так как матрица перемещений – симметричная матрица с убывающими
по обе стороны от диагоналей элементами, выберем 𝛿42 = 𝛿24 =
0,0001, удовлетворяющие полученному выше ограничению.
Результат реализации модели для случая трёх городов в программе
Matlab представлен на Рис.3.
16
Рисунок 3.
Из рисунка видно, что вспышка эпидемии в четвёртом городе
возникает
после
некоторой
временной
задержки.
Количество
инфицированных лиц растёт медленно. Это может быть связано с большой
удалённостью рассматриваемого города от предыдущих и низким уровнем
транспортных перемещений. При этом уровень инфекции остаётся на одном
уровне с предшествующим третьим городом.
17
Выводы
В результате теоретического исследования были найдены условия на
интенсивность транспортных потоков 𝛿𝑖𝑗 , при которых для системы
уравнений, представляющую собой математическую модель эпидемии в
нескольких городах, выполняются условия возникновения эпидемии. Таким
образом,
поставленная
задача
о
нахождении
ограничений
для
𝛿𝑖𝑗 ,
обеспечивающих выполнение условий (2) и (3) была решена полностью.
Решение данной задачи было найдено последовательно для случая двух
городов, а затем для случаев трёх и четырёх. Были заданы конкретные
начальные данные, на основании которых были вычислены ограничения.
Найденные значения интервалов, а также выбранные значения представлены
в Таблице 1.
Коэффициент
Нижнее
Выбранное
Верхнее
пропорциональности
значение 𝛿𝑖𝑗
значение 𝛿𝑖𝑗
значение 𝛿𝑖𝑗
𝛿12 = 𝛿21
0,016
0,003
0,434
𝛿13 = 𝛿31
-0,037
0,002
0,628
𝛿23 = 𝛿32
0
0,001
0,891
𝛿14 = 𝛿41
0,014
0,001
0,245
𝛿24 = 𝛿42
0
0,001
0,508
𝛿34 = 𝛿43
-0,019
0,0001
0,812
Таблица 1.
В
математическом
пакете
Matlab
была
написана
программа,
реализующая движение эпидемической волны для четырёх городов. В
эксперименте
были
использованы
начальные
данные,
выбранные
в
соответствии с требованиями и значения коэффициентов из вычисленных
интервалов.
Численно было показано, что выведенные соотношения для границ
коэффициентов 𝛿𝑖𝑗 соответствуют получаемому в эксперименте.
18
Заключение
Построенная
математическая
модель
описывает
движение
эпидемической волны для случая нескольких городов. Характер и поведение
волн эпидемии во многом зависит от заданных начальных данных. При этом
при изменении начальных условий необходимо производить перерасчёт всех
коэффициентов последовательно для каждого города, начиная с первого.
Основными трудностями, которые возникли при построении данной модели,
являются трудоёмкие вычисления. Было замечено, что объём требуемых
вычислений растёт с увеличением числа городов. Это условие осложняют
применение модели для большого числа городов.
В результате проведённого эксперимента было замечено, что общий
уровень инфицированных лиц в каждом городе уменьшается по мере
удаления от первого города. Кроме того, требуется большее количество
времени для появления вспышки. Данные факты объясняют выбранные
коэффициенты для транспортной матрицы перемещений. То есть, чем
дальше друг от друга расположены города, тем меньше они связаны между
собой транспортными потоками. В соответствие с чем уровень передаваемой
инфекции более низок.
19
Список литературы
1. Смородинцев А.А. Грипп и его профилактика (Руководство для
врачей). Л.: Медицина, 1984. 384 с.
2. Букринская А.Г. Вирусология. М.: Медицина, 1986. 336 с.
3. Бароян О.В., Рвачев Л.А., Иванников Ю.Г. Моделирование
и
прогнозирование эпидемий гриппа для территории СССР. М.:ИЭМ
им. Н.Ф. Гамалеи, 1977. 546 с.
4. Kermack W.O., McKendrick A.G. A Contribution to the Mathematical
Theory of Epidemics // Proceedings of the Royal Society of London. Series
A, Containing Papers of a Mathematical and Physical Character, 1927. Vol.
115, No 772. P. 700-721.
5. Колесин И.Д., Житкова Е.М. Математические модели эпидемий.
Учебное пособие. СПб: НИИФ СПбГУ, 2004. 92 с.
20
Приложение
Ниже
представлен
код
программы
Matlab,
реализующий
математическую модель эпидемии для случая четырёх городов.
function threepoint
clear all; clc;
[T,Y] = ode45(@odefun,[0 100],
[900 100 0 950 50 0 985 15 0 999 1 0]);
figure;
plot (T,Y(:,2),T,Y(:,5),T,Y(:,8),T,Y(:,11),'LineWidth',5);
grid on;
xlabel('t','fontsize',25);
ylabel('I','fontsize',25);
legend ('I1','I2','I3','I4');
function dy=odefun(t,y)
a=0.00025; beta=0.1;
dy=zeros(12,1);
p=[0 300000000 200000000 100000000;
300000000 0 100000000 1000000;
200000000 100000000 0 1000000;
100000000 1000000 1000000 0];
d=0.0000000001*p;
dy(1)=-a*[y(2)+d(2,1)*y(5)+d(3,1)*y(8)+d(4,1)*y(11)]*y(1)d(1,2)*y(1)-d(1,3)*y(1)d(1,4)*y(1)+d(2,1)*y(4)+d(3,1)*y(7)+d(4,1)*y(10);
dy(2)= a*[y(2)+d(2,1)*y(5)+d(3,1)*y(8)+d(4,1)*y(11)]*y(1)d(1,2)*y(2)-d(1,3)*y(2)d(1,4)*y(2)+d(2,1)*y(5)+d(3,1)*y(8)+d(4,1)*y(11)-beta*y(2);
dy(3)= beta*y(2)-d(1,2)*y(3)-d(1,3)*y(3)d(1,4)*y(3)+d(2,1)*y(6)+d(3,1)*y(9)+d(4,1)*y(12);
dy(4)=-a*[y(5)+d(1,2)*y(2)+d(3,2)*y(8)+d(4,2)*y(11)]*y(4)d(2,1)*y(4)-d(2,3)*y(4)d(2,4)*y(4)+d(1,2)*y(1)+d(3,2)*y(7)+d(4,2)*y(10);
dy(5)= a*[y(5)+d(1,2)*y(2)+d(3,2)*y(8)+d(4,2)*y(11)]*y(4)d(2,1)*y(5)-d(2,3)*y(5)d(2,4)*y(5)+d(1,2)*y(2)+d(3,2)*y(8)+d(4,2)*y(11)-beta*y(5);
dy(6)= beta*y(5)-d(2,1)*y(6)-d(2,3)*y(6)d(2,4)*y(6)+d(1,2)*y(3)+d(3,2)*y(9)+d(4,2)*y(12);
dy(7)=-a*[y(8)+d(1,3)*y(2)+d(2,3)*y(5)+d(4,3)*y(11)]*y(7)d(3,1)*y(7)-d(3,2)*y(7)d(3,4)*y(7)+d(1,3)*y(1)+d(2,3)*y(4)+d(4,3)*y(10);
dy(8)= a*[y(8)+d(1,3)*y(2)+d(2,3)*y(5)+d(4,3)*y(11)]*y(7)d(3,1)*y(8)-d(3,2)*y(8)d(3,4)*y(8)+d(1,3)*y(2)+d(2,3)*y(5)+d(4,3)*y(11)-beta*y(8);
dy(9)= beta*y(8)-d(3,1)*y(9)-d(3,2)*y(9)d(3,4)*y(9)+d(1,3)*y(3)+d(2,3)*y(6)+d(4,3)*y(12);
21
dy(10)=-a*[y(11)+d(1,4)*y(2)+d(2,4)*y(5)+d(3,4)*y(8)]*y(10)d(4,1)*y(10)-d(4,2)*y(10)d(4,3)*y(10)+d(1,4)*y(1)+d(2,4)*y(4)+d(3,4)*y(7);
dy(11)= a*[y(11)+d(1,4)*y(2)+d(2,4)*y(5)+d(3,4)*y(8)]*y(10)d(4,1)*y(11)-d(4,2)*y(11)d(4,3)*y(11)+d(1,4)*y(2)+d(2,4)*y(5)+d(3,4)*y(8)-beta*y(11);
dy(12)= beta*y(11)-d(4,1)*y(12)-d(4,2)*y(12)d(4,3)*y(12)+d(1,4)*y(3)+d(2,4)*y(6)+d(3,4)*y(9);
end
end
22
Отзывы:
Авторизуйтесь, чтобы оставить отзыв