CC BY
20
КРИЗИС 2020 г. - ПАНДЕМИЯ КОРОНАВИРУСА
Д.М. Ксенофонтов
СЦЕНАРНОЕ МОДЕЛИРОВАНИЕ ЭПИДЕМИОЛОГИЧЕСКИХ ЭФФЕКТОВ ЭКОНОМИЧЕСКОЙ ПОЛИТИКИ
Б01: 10.47711/2076-318-2020-542-565
Введение. Пандемия С0УГО-2019 показала необходимость разработки точных моделей развития эпидемии, которые могли бы ответить на вопрос о том, когда и какие меры должны быть приняты государствами для предотвращения катастрофы [1]. Однако очевидно, что решение о введении карантина, останавливающего целые сектора экономики, но при этом спасающего жизни -это только одна сторона медали. Карантин не может длиться неограниченно долго, рано или поздно необходимо принимать решение о возвращении к нормальной жизни, так как негативные социально-экономические последствия остановки экономики могут перевесить медицинские опасности, связанные с распространением инфекции.
Для принятия взвешенных решений о продлении или снятии карантинных мероприятий необходим всесторонний анализ как экономических последствий продления, так и медицинских эффектов снятия ограничений, так как очевидно, что запуск остановленных предприятий и возвращение к обычной жизни неизбежно приведет к росту числа контактов и потенциально к разгону эпидемии. Главной задачей прогнозов, предоставляемых медицинскими моделями в этом случае, является оценка сценариев запуска экономики с точки зрения возможной перегрузки системы здравоохранения. Таким образом, модели развития эпидемии должны стать неотъемлемой частью обоснования решений экономической политики в условиях пандемии.
Целью настоящего исследования является отработка методики взаимодействия экономических и медицинских моделей для совместного анализа решений экономической политики в условиях эпидемии С0УГО-2019. Достижение поставленной цели требует разработки математической модели эпидемии и инструментов аппрокси-
542
мации, которые могли бы адекватно описать ретроспективные данные о динамике эпидемии в отдельном регионе, оставаясь при этом сравнительно вычислительно простыми (не требующими суперкомпьютерных мощностей), а также отработки взаимодействия и обмена данными с экономическим блоком моделей.
Математическая модель эпидемии. Главная трудность определения параметров модели эпидемии СОУГО-19 по статистическим данным связана с принципиальным различием между числом выявленных случаев и общем числом активных зараженных. При этом нужно учитывать, что выявленные зараженные с высокой вероятностью изолируются и не будут участвовать в дальнейшем распространении инфекции, что существенно влияет на скорость распространения инфекции в популяции. Поэтому, модель, способная корректно аппроксимировать динамику распространения заболевания, отраженную в официальной статистике, должна состоять из двух частей - скрытой и наблюдаемой, связанных между собой через выявляемость, возможно нестационарную.
Математическая модель, разработанная в рамках данного исследования, должна была учесть ряд факторов: во-первых, у инфекции есть инкубационный период, который вносит задержку между изменением условий распространения и отражением этих изменений в статистике. Далее, модель должна отражать существование пресимптомной стадии течения заболевания, когда инфицированный уже заразен, но явных симптомов еще нет [2; 3], так как наличие таких «супер-распространителей» может существенно сказаться на динамике эпидемии. Также, модель должна различать легкие (не угрожающие жизни и допускающие амбулаторное лечение) и тяжелые (требующие госпитализации) формы течения заболевания, с возможностью перехода зараженного с тяжелой формой в критическое состояние, требующего интенсивной терапии, что позволяет учесть ограниченность ресурсов медицинской системы. Наконец, модель должна учитывать интенсивность выявления и изоляции (выявляемость) зараженных с различными формами течения заболевания для обеспечения связи между скрытой и наблюдаемой частями. Стоит отметить, что выявляемость может учитывать не только медицинские факторы, такие как объем тестирования, но и психологические. К примеру, на фоне страха перед эпидемией, человек с симптомами ОРВИ с большей вероятностью обратится в медицинское учреждение или изолируется.
543
Для сценарного моделирования развития эпидемии была построена модель среднего поля, являющаяся развитием классической модели SEIR [4]. В рамках разработанной модели популяция региона была разделена на следующие субпопуляции (компартменты, объединение всех субъектов, находящихся в одном состоянии):
• восприимчивые к инфицированию (5);
• контактные (E);
• латентные / пресимптомные (IL, ILo);
/ tM TM \
• легкие (1 );
• тяжелые (I5);
• критически тяжелые (Iе);
• погибшие (D);
• выздоровевшие (R).
Индексом ISO обозначены выявленные и изолированные зараженные.
Модель получила наименование SEIRD-LMSC-ISO по моделируемым субпопуляциям, как это принято для моделей среднего поля, которое расшифровывается как «Восприимчивые(5)-Контактные(£)-Зараженные(Т)-Выздоровевшие^-Пошбшиеф) с разделением зараженных на пресимптомных^), легких(М), тяжелых(5) и критически тяжелых (C) и учетом выявления и изоляции(КО)».
Схема переходов между субпопуляциями представлена на рис. 1. Каждому переходу соответствует определенная интенсивность (число событий в единицу времени).
Рис. 1. Схема переходов модели SEIRD-LMSC-ISO
544
Предполагается, что зараженные с тяжелой формой являются выявленными и изолированными и все летальные исходы связаны с заболеванием.
К наблюдаемой части модели относятся следующие показатели: суммарное число выявленных пресимптомных ), легких форм
(1Ц0) и тяжелых форм (I5), а так же число летальных исходов (Б).
Число выздоровевших к наблюдаемым показателям не относится, так как модель и статистика трактуют этот показатель по-разному (фактическое прекращение периода заразности и выписка из медицинского учреждения/прекращение медицинского наблюдения соответственно).
Эволюция численностей субпопуляций описывается системой нелинейных дифференциальных уравнений (1)-(18)\ Описание
параметров системы представлено в табл. 1.
5 = -в(г)(11гее /Ы)■ 5 (1)
Е = в(г)(1^е / N) ■ 5 -5Е^ ■ Е (2)
Iь = 5е^ ■ Е -5^ ■ I, -0ь(г) ■ I, (3)
¡г =0ь(г)■ 1Ь -5^ ■ II0 (4)
1 м = Р^ы5ь^лЬ - 0ы(*) ■ ^ - Ум ■ ^ (5)
^ = Рь^м 5 ^Л0 + 0м (() ■ ^ - Ум ■ РМ0 (6)
^ = (1 -Рь^м)5ь^л ■ IЬ -(1 -Р^Л5 ■ ^ -Р^с5^ ■ ^ (7)
I с = Р^с55^с ■ I, - (1 -а)1 с ■ min(Cmed'Ic) - (8) - ар ■ тт(стеЛ,Ic) -тах(0,^ - стеЛ )
Б = ар ■ тт( С^,^ ) + тах( 0(9)
Р = Ум ■ II + (1 -Р5^с)У5 ■ ^ + (1 - а)ус ■ ттСт^.Ь) (10)
I* = I, + I = Ь,м (11)
I*" = ^ + Км^м + X ■ (Iь + к); К, х < 1 (12)
Р(0 = уе# -ад (13)
ъ =1/Г^к, I = м, 5, с (14)
(15)
1 Точка над переменной означает дифференцирование по времени.
545
^L^-A 1/Tactivation (16)
=1/Tcritical (17)
P= 1/Tdeath (18)
Таблица 1
Параметры модели SEIRD-LMSC-ISO
№ Параметр Описание Класс параметра
1 Ro(t) Базовое репродуктивное число. Данный параметр равен числу человек, которые заражаются одним инфицированным до его изоляции или прекращения заразности. Скорость распространения
2 TM^R Временной период сохранения заразности инфицированного с легкой формой течения заболевания Скорость распространения
3 Ts^R Временной период выздоровления зараженного с тяжелой формой течения заболевания Загруженность медицинской системы
4 TC^R Временной период выздоровления зараженного в критически тяжелом состоянии Загруженность медицинской системы
5 Tincubation Инкубационный период - период развития заразности. Скорость распространения
6 Tactivation Пресимптомный период - временной интервал развития заметных симптомов, при этом инфицированный уже заразен. Скорость распространения
7 Tcritical Временной период развития критически тяжелых симптомов у зараженного с тяжелой формой течения заболевания Загруженность медицинской системы и динамика смертности
8 W Выявляемость зараженных с пресимптомной стадией заболевания. Скорость распространения
9 6m(0 Выявляемость зараженных с легкой формой заболевания Скорость распространения
10 pL^M Вероятность легкой формы течения заболевания Скорость распространения
11 PS^C Вероятность критически тяжелой формы течения заболевания Загруженность медицинской системы
12 a Вероятность летального исхода у зараженного в критически тяжелом состоянии. Динамика смертности
13 Tdeath Временной период от появления критически тяжелых симптомов до летального исхода Динамика смертности
14 Kl Доля выявленных пресимптомных зараженных, продолжающих распространять инфекцию Скорость распространения
15 K m Доля выявленных зараженных с легкой формой, продолжающих распространять инфекцию Скорость распространения
16 K Понижающий коэффициент заразности пресимптомных инфицированных по сравнению с симптомными Скорость распространения
17 Cmed Емкость медицинской системы - число пациентов в критически тяжелом состоянии, которому одновременно может оказываться помощь. Загруженность медицинской системы и динамика смертности
546
Каждый параметр данной модели может быть задан функцией от времени и численностей субпопуляций, что позволяет исследовать не только инерционные сценарии, но и влияние принимаемых мер, а также получать самосогласованные решения, такие как адаптивный карантин.
Разработанная модель имеет 17 свободных параметров, если коэффициенты считаются постоянными во времени и независимыми от численностей субпопуляций.
Модель была реализована в виде пакета программ на языке Python 3. Интегрирование дифференциальных уравнений выполнялось методом Рунге-Кутты с помощью библиотеки ODEPACK, для обеспечения высокой производительности. Шаг выдачи результатов задавался равным одному дню, а шаг интегрирования подбирался библиотекой автоматически для обеспечения устойчивости явной схемы интегрирования по времени.
Данная программная реализация позволяет гибко перенастраивать модель, передавая пользовательские функции заданной сигнатуры (заданного числа позиционных аргументов) для вычисления любых коэффициентов правой части системы уравнений.
В рамках этого исследования субпопуляции не делились на возрастные группы несмотря на то, что такая векторизация поддерживается программной реализацией, так как это существенно увеличивает количество параметров для аппроксимации. Параметры эпидемии принимались усредненными по всем возрастным группам.
Оценка параметров модели SEIRD-LMSC-ISO на основе аппроксимации ретроспективных данных о распространении COVID-2019 в Москве. В силу принципиальной неоднородности параметров распространения и стадии эпидемии в регионах России, для анализа была выбрана Москва, где из-за высокой плотности населения и большой скорости распространения условия наиболее близки к допущениям моделей среднего поля - «термодинамичность» и однородность моделируемой популяции. Термодинамичность системы предполагает, что численность и плотность населения достаточно высока, чтобы можно было использовать центральную предельную теорему и перейти от вероятностного описания процесса распространения инфекции к непрерывному и детерминированному. Однородность популяции означает, что каждый субъект имеет одинаковую вероятность заражения. Это предположение работает в тех случаях, когда присутст-
547
вует интенсивное перемещение субъектов с большим числом контактов. Для Москвы это допущение обеспечивается наличием развитой системы общественного транспорта, включая метро.
Для аппроксимации были выбраны ретроспективные данные о суммарном числе выявленных случаев и полном числе летальных исходов в Москве с 06.03.2020 по 11.05.2020. Данные с 12.05.2020 по 22.05.2020 использовались в качестве контрольных, по которым оценивалось качество краткосрочного прогноза.
Для учета мер противодействия распространению инфекции, R0, eL и Ум задавались как непрерывно-дифференцируемые функции от времени. Зависимость базового репродуктивного числа (19) от времени состояла из трех условных интервалов: распространение без противодействия (t< tR), адаптация популяции к карантину (tR < t< tR)
и устойчивый карантин (t> t1R). Каждому из интервалов соответствовало свое значение базового репродуктивного числа - R0tart, Rm'ddle и R0"d соответственно. Скорость перехода от одного значения к другому определялось параметром kR большие значения которого соответствуют более быстрому переходу.
R0 = (ЯГ' -R™ddk)/{ 1 + exp[-kR • (-t + tR)]} + + (Rm,Me -Rn)/{1 + exp[-kR .(-t + tR)]} + Rend
Зависимость выявляемости от времени задавалась выражением (20) и состояла из двух условных интервалов - первому соот-
0 start end
M , второму - УM , скорость перехода между значениями интервалов определяется параметром ke. При проведении аппроксимации использовалась только выявляемость легких форм, выявляемость же зараженных на пресимптомной стадии вычислялась как Ум, умноженная на понижающий коэффициент (0.5).
Ум = (УГ-ee0nd)/{1 + exp[-ky (-t + ty)]} + eeMd. (20)
Помимо параметров модели, необходимо было определить адекватные начальные условия, которые обеспечивают скорость распространения, хорошо совпадающую со статистикой и нахождение базового репродуктивного числа в адекватном диапазоне (2-7) [5]. В качестве начального условия для каждого расчета задавалось определенное количество зараженных на пресимптомной стадии, которое также было варьируемым параметром.
548
С учетом временных зависимостей (19)-(20) и необходимости определить начальное условие, для аппроксимации ретроспективных данных требуется определение 23 свободных параметров. Такая конфигурация пространства параметров делает невозможным поиск решения методом «грубой силы» - прямого перебора параметров с заданным шагом. Более того, прямым перебором невозможен даже поиск качественного начального приближения для оптимизации.
Для получения решения был реализован метод случайного поиска начального приближения с последующей минимизацией невязки методом последовательного квадратичного программирования (SLSQP/RSS) [6]. Одна итерация цикла решения состояла из двух шагов: на первом шаге из заданного для каждого параметра диапазона выбиралось случайное начальное значение, на втором - выполнялась минимизация невязки из заданной на первом шаге точки пространства параметров. Невязка оценивалась в квадратичной норме по нормированным на максимальное значение ретроспективным кривым динамики полной смертности и суммарного числа выявленных зараженных для Москвы. Вес невязки по каждой из кривых равнялся единице.
Программная реализация была выполнена на языке Python 3 с применением библиотек numpy/scipy и модуля multiprocessing.
Результат каждой итерации цикла решения сохранялся для последующей классификации и анализа. Всего было проведено 16000 итераций цикла решения, что потребовало 14 часов вычислений на 24-процессорных ядрах.
После окончания поиска решения, была отобрана первая квартиль по величине невязок. На основе отобранных решений были построены доверительные интервалы (с вероятностью 99,7%) и средняя кривая. Затем, были выбраны параметры, дающие решение с наименьшей невязкой. Оптимальные параметры, полученные в результате обработки результатов представлены в табл. 2.
Качество аппроксимации ретроспективных данных о ежедневном приросте выявленных случаев и полном числе летальных исходов, связанных с заболеванием можно оценить на рис. 2, 3.
Для оптимальной модели, коэффициент детерминации ежедневного прироста выявленных случаев, рассчитанный по аппроксимируемым данным, составил 0,987.
549
Таблица 2
Диапазоны поиска и параметры наилучшей аппроксимации
№ Параметр Диапазон изменения Значение наилучшей аппроксимации
1 IL 10 100-10000 240
2 restart R0 2-7 6.7
3 Т") middle R0 0.5-3.5 0.55
4 г* end R0 0.5-3.5 1.83
5 kR 0.1-5.0 0.19
6 tR l0 4.04.2020-26.04.2020 10.04.2020
7 tR 1 27.04.2020-7.05.2020 5.05.2020
8 TM^R 5.8-10.8 8.7
9 Ts^R 6.9-12.9 7.0
10 Tc^R 6.7-12.5 7.7
11 Tincubation 3.8-7.0 4.2
12 Tactivation 2.3-4.3 4.1
13 Tcritical 6.0-12.0 11.4
14 r\ start e M 0.001-0.1 0.0076
15 /л end e M 0.01-0.3 0.26
16 ke 0.05-5.0 0.31
17 19 17.04.2020-7.05.2020 5.05.2020
18 PL^M 0.59-0.81 0.79
19 Ps^C 0.1-0.8 0.1
20 a 0.4-0.8 0.63
21 Tdeath 10.6-19.8 13.0
22 KL 0.015-0.085 0.062
23 Km 0.0075-0.0425 0.021
550
Рис. 2. Ежедневный прирост выявленных случаев:
-Модель (средняя кривая); ♦ данные для аппроксимации;
О контрольные данные; □ доверительный интервал 99,7%
4500 4000 -3500 3000 2500 -2000 1500 1000 -500 0
Период
Рис. 3. Полное число летальных исходов: -Модель (средняя кривая); ♦ данные для аппроксимации;
О контрольные данные; □ доверительный интервал 99,7%
Аналогичный параметр для суммарного числа летальных исходов равнялся 0,992. Контрольные данные для обеих кривых укладываются в доверительный интервал 99,7%.
Как можно видеть, контрольные данные для полного числа летальных исходов смещаются к верхней границе доверитель-
551
ного интервала, т.е. модель имеет тенденцию к занижению смертности, при том что контрольные данные по ежедневному приросту выявленных зараженных располагаются близко к средней кривой (средняя невязка не превышает 12%). Такое поведение связано с принципиальной разницей в аппроксимации производной (ежедневный прирост числа заболевших) и интегральной величины - числа летальных исходов. Очевидно, что в случае аппроксимации интегрального показателя, точность вычисления производной на границе отрезка будет существенно ниже, если не применять дополнительные меры, такие как введение нестационарных весов для невязки, возрастающих к границе области аппроксимации. Выбор адекватной взвешивающей функции является самостоятельной задачей и выходит за границы данного исследования.
Рассмотрим, как соотносятся динамика числа активных случаев в скрытой и наблюдаемой частях модели для оптимальных параметров (рис. 4).
Число активных случаев
100000 90000 80000 70000 60000 50000 40000 30000 20000 10000 0
с^т^^'ЛЮЮ^сосо Период
оооооооооо
сКосКсКсКсососо^!^
Рис. 4. Динамика числа активных случаев в наблюдаемой части модели (-)
и полное число активных случаев (-) для оптимальных параметров
Можно видеть, что наблюдаемый пик эпидемии, связанный с ростом выявляемости, отстает от истинного.
Стоит подчеркнуть, что, несмотря на достаточно хорошую аппроксимацию ретроспективных данных, данная модель является инструментом сценарного моделирования, а не универсальным средством прогнозирования развития эпидемии.
552
Исследование особенностей модели БЕТКЭ-ЬМБС-ШО. Рассмотрим ряд сценариев, которые наиболее полно раскрывают особенности и возможности разработанной модели с использованием оптимальных параметров.
С точки зрения сценарного моделирования экономических последствий эпидемии, или, наоборот, медицинских эффектов социально-экономических решений, наиболее важными являются такие «низкочастотные» показатели, как, например, смертность в годовом исчислении и нагрузка на медицинскую систему.
Инерционный сценарий развития эпидемии.
Для начала, рассмотрим инерционный сценарий, когда меры противодействия распространению не применяются. В рамках данного сценария, базовое репродуктивное число было принято постоянным во времени и равнялось 5. Также предполагалось, что выявляются только тяжелые формы заболевания. Остальные параметры равнялись представленным в табл. 2. Начальное условие - 1000 контактных. На рис. 5 представлен результат моделирования смертности (отношение числа летальных исходов к населению), связанной с С0УГО-2019 в Москве в годовом исчислении для различных емкостей медицинской системы, которая в рамках данной серии расчетов представляет собой максимальное число зараженных в критически тяжелом состоянии, которым может одновременно оказываться помощь, безотносительно к числу аппаратов ИВЛ и другого специального оборудования.
Смертность
3,0% 2,5% 2,0% 1,5% 1,0% 0,5% 0,0%
1,6%
неограниченная 45000 30000
Количество коек в 15000 больницах
2,5%
1.9%
1,1%
Рис. 5. Смертность в годовом исчислении, связанная с С0У1Б-2019 в Москве при отсутствии мер противодействия
553
Как можно видеть, даже в случае с неограниченной емкостью медицинской системы, отсутствие мер противодействия приводит к смертности, связанной с заболеванием, сравнимой со смертностью от всех причин в Москве в 2019 г (0,95%, ~ 120000 человек). При наличии ограничений на емкость больниц, смертность от С0УГО-2019 может возрастать в 2,5 раза и более (при уменьшении емкости системы здравоохранения). Стоит учитывать, что модель имеет тенденцию к занижению смертности, так что показатели инерционного сценария могут быть еще хуже. Конечно, инерционный сценарий кажется слабо реализуемым на практике, так как даже если меры не принимаются на уровне государства, население, на фоне страха перед эпидемией, самостоятельно ограничит контакты и, скорее всего, замедлит распространение инфекции, что можно наблюдать, например, в Швеции. Тем не менее, данный сценарий однозначно показывает необходимость противодействия распространению такого высоко контагиозной и сравнительно летальной инфекции, как 8ЛЯ8-иС0У-2. Влияние задержки объявления карантина на характер течения эпидемии
Оценим влияние задержки объявления карантина на пиковую загрузку больниц и смертность в годовом исчислении. Пусть емкость системы здравоохранения составляет 15000 коек. В рамках данной серии расчетов, базовое репродуктивное число до объявления карантина равнялось 5, во время карантина - 0,95. Предполагалось, что выявляются только тяжелые формы заболевания. Эпидемия начиналась с 1000 контактных.
Рассмотрим, как влияет задержка объявления карантина на пиковую нагрузку на систему здравоохранения (рис. 6). Каждое значение на рис. 6 представляет собой максимальную загрузку больниц, достигнутую в рамках отдельного сценария развития эпидемии, при котором карантин был объявлен с задержкой, отложенной по оси абсцисс.
Данный график показывает, что существует определенная задержка объявления карантина, после которой перегрузка системы здравоохранения неизбежно наступает, что должно неизбежно привести к росту, как ежедневного прироста летальных исходов, так и смертности в годовом исчислении. Для базового репродуктивного числа, до карантина равного 5 и 1000 контактных, эта критическая задержка составляет 55 дней.
554
Пиковая загрузка коечного фонда, %
Рис. 6. Зависимость пиковой нагрузки на систему здравоохранения (загрузки коечного фонда) от задержки объявления карантина
Рассмотрим, как длительность карантина влияет на пиковый ежедневный прирост летальных исходов и смертность, связанную с СОУГО-2019, в годовом исчислении (рис. 7).
Пиковая дневная Смертность, %
смертность
Рис. 7. Зависимость смертности (---), связанной с заболеванием,
в годовом исчислении и пикового дневного прироста летальных исходов
(пиковой дневной смертности -) от задержки объявления карантина
(левая ось ординат имеет логарифмический масштаб)
Необходимо отметить, что до перегрузки системы здравоохранения смертность (как в дневном, так и в годовом исчислении) демонстрируют экспоненциальный рост с задержкой. Таким образом, с медицинской точки зрения, скорость реакции играет определяющую роль на демографический эффект эпидемии. При
555
достижении критической задержки объявления карантина в 55 дней, при которой гарантированно наступает перегрузка системы здравоохранения, смертность начинает быстрее расти с задержкой. При дальнейшем увеличении задержки, смертность стремится к асимптотическому значению, равному смертности при отсутствии мер противодействия распространению инфекции и карантинные меры становятся бессмысленными. Это связано с тем, что пики численности всех субпопуляций уже пройдены к этому моменту и самый пессимистичный сценарий эпидемии уже реализован. Влияние длительности карантина и выявляемости на характер течения эпидемии. Эффекты разгона эпидемии. Оценим влияние длительности карантина на смертность, связанную с заболеванием и нагрузку на систему здравоохранения. В рамках данной серии расчетов емкость больниц равнялась 15000 койкам, эпидемия начиналась с 1000 контактных, карантин объявлялся на 40 день эпидемии, а также считалось, что выявляются только зараженные с тяжелыми формами. Параметр Я0 до карантина равнялся 5, во время - 0,95, после - 1,85. Остальные параметры модели совпадали с представленными в табл. 2.
Как видно из рис. 8-9, чувствительность смертности и нагрузки на систему здравоохранения к изменению длительности карантина существенно ниже чувствительности к задержке объявления карантина - при изменении длительности карантина фактически на 2 порядка, характеристики эпидемии варьируются в диапазоне от 20 до 30%.
Пиковая дневная
смертность
Смертность, %
750
0,70
700
650
0,65
600
0,60
550
500
0,55
Длительность карантина, дни
10
30
50
70
90
Рис. 8. Зависимость смертности (---) в годовом исчислении
и пиковой дневной смертности (-) от длительности карантина
556
Пиковая загрузка коечного фонда, %
100 -
95 -90 -85 -80 -
10
30
50
70
90
Длительность карантина, дни
Рис. 9. Зависимость пиковой нагрузки на систему здравоохранения от длительности карантина
Это легко объяснимо, если рассмотреть динамику критически тяжелых случаев (рис. 10) для нескольких длительностей карантина (30, 60 и 90 дней).
Число критически тяжелых
14000 12000 -10000 -8000 -6000 4000 2000 0
Дни
Рис. 10. Зависимость числа критически тяжелых случаев от времени для различных длительностей карантина: -30 дней;---60 дней; -----90 дней
Как можно видеть, основная часть критически тяжелых случаев и, следовательно, летальных исходов связана с разгоном эпидемии после снятия карантина. Это одна из особенностей математической формулировки классических моделей среднего поля - если базовое репродуктивное число больше 1, а в по-
557
пуляции есть восприимчивые и зараженные, эпидемия будет развиваться. Рассмотрим, как изменится динамика эпидемии, если модель учитывает выявляемость и изоляцию легких форм течения СОУГО-19. Предположим, что на сотый день эпидемии начинается масштабное тестирование населения с выявляемо-стью зараженных с легкой формой заболевания 9% в день. Как можно видеть (рис. 11), наличие выявляемости зараженных с легкой формой течения заболевания радикально изменяет картину развития эпидемии - даже для длительности карантина в 30 дней пик связанный с возвращением к обычной жизни меньше в 3.5 раза меньше, чем в сценарии без выявляемости. А для длительности 60 дней и разгона эпидемии после снятия карантина вообще не происходит. Таким образом, выяв-ляемость (и изоляция выявленных) зараженных играет сущее-ственную роль в сдерживании эпидемии, наравне с карантинными мероприятиями.
Число критически тяжелых
14000 -12000 -10000 -8000 -6000 -4000 2000 0 -■
Рис. 11. Зависимость числа критически тяжелых случаев от времени для различных длительностей карантина при наличии выявляемости зараженных с легкой формой течения заболевания Интенсивность выявления - 9% в день -30 дней;---60 дней (Я0=0,95); -— 90 дней (Л0=0,95)
Взаимодействие с блоком расчета экономических последствий
Экономические последствия карантинных мероприятий рассчитывались с помощью модели-мультипликатора, управляемыми параметрами которой являются даты остановки и запуска отраслей эко-
558
номики моделируемого региона . Предполагалось, что 33 из 59-ти отраслей были остановлены 30.03.2020 г. Остальные 26 отраслей продолжили работать. Даты запуска остановленных отраслей варьировались в диапазоне от 04.05 до 15.07 в зависимости от отрасли. Прямой эффект сокращения ВВП в данном сценарии составил 7,3%.
Влияние карантинных мер на отрасли экономики в данном сценарии представлено на рис. 12-13.
0% -2% -4% -6% -8% 0% 2% 4% 6%
И
Код ОКВЭД 1 Л
Рис. 12. Виды деятельности со спадом производства более 6%
1% , 0% --1% -2% --3% -4% -5% --6%
^'л'люо^^г^'л^гс^^ю'лтс^то^ос^с^'л^гоото^от^ Код ОКВЭД 1.1
Рис. 13. Виды деятельности со спадом производства менее 6%
Выражаю глубокую благодарность д.э.н. Александру Александровичу Широву за предоставленную для отработки методологии модель экономических эффектов карантина.
559
Схема взаимодействия расчетных блоков представлена на рис. 14. Основным связующим параметром является прирост численности работающего населения, обусловленный запуском отраслей экономики. Очевидно, что возрастание числа экономически и социально активных субъектов в популяции приводит к неизбежному росту числа контактов и увеличению репродуктивного числа. Рост Я0 с изменением числа свободных субъектов рассчитывается медицинским блоком.
Рис. 14. Схема взаимодействия экономического и медицинского расчетных блоков
560
Оценка медицинских эффектов сценариев запуска отраслей экономики Москвы.
Оценим медицинские эффекты запуска отраслей экономики на основе базового экономического сценария снятия ограничений и оптимальный параметров модели распространения COVID-2019 в Москве. Входными данными для SEIRD-LMSC-ISO в данной постановке задачи является доля трудоспособного населения, вернувшегося к работе (рис. 14). Изменение базового репродуктивного числа, связанного с выходом из карантина, носила аддитивный характер [7].
Связь между экономической и медицинской частью системы обеспечивалась с помощью библиотеки xlrd, позволяющей получать доступ к значениям ячеек в файлах xlsx/xlsm из программ, написанных на языке Python 3. Динамика выхода работающего населения из карантина, полученная от экономической части системы, представлена на рис. 15.
Доля работающего населения
Рис. 15. Динамика выхода работающего населения из карантина, соответствующая базовому сценарию запуска отраслей
На рис. 16-17 представлено сравнение динамики ежедневного прироста выявленных случаев и полного числа летальных исходов соответственно для сценариев снятия и сохранения карантина.
561
7000 6000 -5000 -4000 -3000 -2000 -1000 -0
Сценарии:
- снятия карантина
- сохранения карантина
Период
Рис. 16. Динамика ежедневного прироста выявленных зараженных для сценариев сохранения и снятия карантина: • данные для аппроксимации; О контрольные данные
4000 3500 3000 2500 2000 1500 1000 500 0
Сценарии: снятия карантина
сохранения карантина
Период
Рис. 17. Динамика полного числа летальных исходов
для сценариев сохранения и снятия карантина: • данные для аппроксимации; О контрольные данные
Как можно видеть, модель предсказывает появление медленно убывающей «полки» прироста новых случаев при снятии карантина, но сохранении выявляемости на уровне оптимальных параметров аппроксимации.
Для оценки медицинских последствий снятия карантина более показательны чистые эффекты - приросты показателей по срав-
562
нению со сценарием сохранения карантина. Снятие карантина в рамках данной модели приводит к тому, что дополнительно заразятся 94,5 тыс. чел. (на 27% больше по сравнению со сценарием сохранения карантина), при этом добавленный максимальный прирост подтвержденных случаев составит 13,4% значения наблюдаемого пика эпидемии (860 чел./сут.). Максимальная добавленная нагрузка на отделения интенсивной терапии составит 15% от пиковой (150 пациентов), а полное число летальных исходов возрастет на 21% по сравнению со сценарием сохранения карантина (700 летальных исходов). В рамках данного сценария снятия карантина, система здравоохранения далека от перегрузки, так как число критически тяжелых пациентов на момент начала снятия карантина составляет 600 и падает до 200 к концу прогнозного периода, что значительно меньше числа коек интенсивной терапии, оснащенных аппаратами искусственной вентиляции легких (около 5000 для Москвы, согласно данным Минздрава России).
Заключение. В рамках данного исследования были отработаны методологические и технические аспекты построения экономико-медицинской системы моделей для сценарного прогнозирования эффектов карантина. Основой медицинского расчетного блока стала разработанная модель среднего поля 8Е1К0-ЬМ8С-180, которая позволила учесть скрытую и наблюдаемую динамику эпидемического процесса, выявляемость различных форм течения заболевания и ограниченную емкость системы здравоохранения. При этом созданная модель отличается невысокой вычислительной сложностью (по сравнению с агентными стохастическими моделями), что позволило, несмотря на большое количество неизвестных, построить алгоритм исследования пространства параметров для аппроксимации ретроспективных данных. Разработанные инструменты аппроксимации могут представлять самостоятельный интерес, так как позволяют не только определить параметры развития эпидемии на ретроспективе, но и построить статистически обоснованную область возможных траекторий развития ситуации при постоянных параметрах модели. Для отработки методики аппроксимации использовались официальные данные о распространения коронавируса в Москве с 6 марта по 12 мая 2020 г. Контроль качества аппроксимации был выполнен на данных с 13 по 28 мая, при этом контрольные данные уложились в доверительный интервал, построенный для вероятности 99,7%.
563
На основе оптимальных параметров модели был исследован сценарий распространения С0УГО-2019 в Москве без мер противодействия (инерционный сценарий). Результаты расчетов демонстрируют смертность (в годовом исчислении), связанную с заболеванием, достигающую 2,5%, что обосновывает необходимость карантинных мероприятий.
Исследовано влияние задержки объявления карантина на характеристики эпидемии (нагрузка на систему здравоохранения и смертность, связанную с заболеванием) и показано, что скорость реакции имеет определяющее значение, так как существует критическая задержка объявления карантина, при которой неизбежно наступает перегрузка медицинской системы, приводящая к резкому росту смертности и в пределе стремящаяся к смертности инерционного сценария.
Оценено влияние выявляемости зараженных с легкой формой течения заболевания на динамику эпидемии при снятии карантина. Показано, что без мер выявления зараженных, эпидемия начинает разгоняться, если к моменту окончания карантина в популяции сохранились активные зараженные. Однако, выявляемость легких форм позволяет стабилизировать динамику заражений и избежать второго пика эпидемии, связанного с возвращением к нормальной жизни. Таким образом, меры по раннему выявлению и изоляции зараженных являются эффективным инструментом противодействия распространению инфекции, наравне с карантинными мероприятиями.
Отработано взаимодействие экономического и медицинского расчетных блоков на примере моделирования медицинских эффектов снятия карантина в соответствии со сценарием остановки / запуска отраслей экономики, предполагающим спад экономики в 7%. Связь моделей осуществлялась через зависимость базового репродуктивного числа от числа трудоспособного населения Москвы, вышедшего из карантина. Расчеты показали, что в рамках разработанной математической модели и при использовании параметров, обеспечивающих наилучшее описание ретроспективных данных, сценарий снятия карантина не приводит к разгону эпидемии и значительной дополнительной нагрузке на систему здравоохранения. Согласно результатам моделирования, возвращение к нормальной жизни и запуск отраслей экономики Москвы приводит к максимальной добавленной нагрузке на отделения интенсивной терапии, составляющей 15% пикового значения (150 пациентов), а число критически тяжелых пациентов изменяется от 600 к началу
564
выхода из карантина до 200 к концу прогнозного периода (1 сентября 2020 г.), что существенно ниже числа коек интенсивной терапии, оснащенных аппаратами ИВЛ (около 5000).
Литература и информационные источники
1. Adam D., Special report: The simulations driving the world's response to COVID-19. Nature 580. Рр. 316-318 (2020).
2. Linton NM., Kobayashi T., Yang Y., Hayashi K, Akhmetzhanov A.R., Jung S.-M. Incubation period and other epidemiological characteristics of 2019 novel coronavirus infections with right truncation //A statistical analysis of publicly available case data. J Clin Med. 2020. № 9(2).
3. Li Q., GuanX., Wu P., Wang X., Zhou L., Tong Y. Early transmission dynamics in Wuhan, China, of novel coronavirus-infected pneumonia //N Engl J Med. 2020. № 382(13). Рр. 1199-1207.
4. Kermack, W.O.; McKendrick, A.G. (1927). A Contribution to the Mathematical Theory of Epidemics // Proceedings of the Royal Society A. № 115(772). Рр. 700-721.
5. Linka K., Peirlinck M., Kuhl E., The reproduction number of COVID-19 and its correlation with public health interventions, medRxiv 2020.05.01.20088047
6. Kraft D. A software package for sequential quadratic programming. 1988. Tech. Rep. DFVLR-FB 88-28, DLR German Aerospace Center—Institute for Flight Mechanics, Koln, Germany.
7. Atkeson A. What Will Be the Economic Impact of COVID-19 in the US? Rough Estimates of Disease Scenarios. NBER // Working Paper. № 26867. March 2020.
565
