Mathematical modeling of the spread of COVID-19 in Moscow

Full Text

ВВЕДЕНИЕ В статье для описания распространения эпидемии в Мо- скве используется дискретное логистическое уравнение. Впервые логистическое уравнение в дифференциальной форме применил бельгийский математик Пьер Ферхюльст в 1845 г. [1] для моделирования роста населения. Принципи- альное отличие от математической модели Томаса Мальтуса (представленной в знаменитой работе «Опыт закона о наро- донаселении» [2]), описывающей экспоненциальный рост популяции, заключается в том, что модель Ферхюльста учи- тывала конкуренцию за ресурсы, приводящую к ограничен- ности роста популяции. В 1920 г. логистическое уравнение в дифференциальной форме стало широко использоваться американским биологом Рай- мондом Пирлом [3] для роста народонаселения, который дал этому уравнению фактически вторую жизнь. В настоящее время использование данного уравнения нашло широкое применение в математической биофизике, что хорошо отражено в монографиях российских биофизиков Г.Ю. Ризниченко и А.Б. Рубина [4; 5]. Уравнение Ферхюльста нашло и сейчас применение в работах для моделирования коронавируса COVID-19 в Китае [6], Швеции [7]. Уравнение Ферхюльста применительно к распространению эпиде- мии имеет две неподвижные (стационарные) точки: ȳ1 = 0, ȳ2 = N. Вторая точка является аттрактором (т.е. устойчивой). Стационарной точкой является максимальное значение численности жителей, которые заболеют. C течением вре- dy  dt  y 1  N  y   мени, какое бы значение не принимал показатель роста в уравнении, численность популяции будет стремиться к N только с разным периодом времени выхода на этот стаци- (где параметр λ характеризует скорость размножения; пара- метр N - максимально возможную численность популяции) онар. Существует большая разница между максимальным числом жителей, которые потенциально могут заболеть, и максимальным значением фактически заболевших жите- лей. Первое значение численности много больше численно- сти фактически заболевших жителей. Поэтому в работе используется дискретное логистиче- ское уравнение, которое также имеет два параметра: пока- затель роста численности заболевших, максимальное зна- чение численности жителей, которые потенциально могут быть инфицированы коронавирусом COVID-19. Дискретное логистическое уравнение получило широ- кую известность благодаря работам американского уче- ного М. Фейгенбаума [8], который обнаружил интересные закономерности, получаемые с помощью этого уравнения численность заболевших будет стремиться к нулю, сколь- ко бы заболевших не было в начале. При значениях 1 < α ≤ 3 безразмерная численность попу- ляции заболевших стремится к стационарному устойчивому состоянию х̄, равному [9] image image x   - 1 . (4)  Следовательно, с течением времени численность попу- ляции заболевших в конце эпидемии будет равна   - 1  и создал теорию универсальности для дискретных отобра- жений. В частности, он показал, что при изменении обезраз- y  N   image . (5)   меренного показателя роста в уравнении в интервале 3-4 возникают бифуркации удвоения периода, а при значении показателя роста равным ~3,5699 логистическое уравнение генерирует хаотические колебания. Дискретное логистиче- ское уравнение после этих работ стало широко использо- ваться для моделирования различных процессов [9; 10]. Для адекватного дискретного описания явлений важной характе- ристикой является выбор шага, дискретного интервала вре- мени, на котором рассматриваются численности популяции. По расчетам Нобелевского лауреата Майкла Левитта, сделанного для Китая, за 24 часа 1 человек в Китае мог за- разить 2,2 статистических человека. То есть один человек мог заразиться от другого на интервале времени ~10,9 часов. Поэтому для дискретного логистического уравнения брался интервал времени, равный 12 часам. Этот период близок к 10,9 часам и удобен для сравнения расчетных и фактиче- ских данных. Иными словами численность популяции забо- левших пересчитывается каждые 12 часов, а сверяются дан- ные через сутки. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ РАСПРОСТРАНЕНИЯ КОРОНАВИРУСА COVID-19 Логистическое уравнение для описания распростране- ния коронавируса COVID-19 имеет вид y  y 1 - yn  , (1) Соотношения (4) и (5) верны при сохранении постоянства показателя α на всем протяжении эпидемии. РЕЗУЛЬТАТЫ РАСЧЕТА ПО МАТЕМАТИЧЕСКОЙ МОДЕЛИ Для проверки адекватности использования логистиче- ского уравнения распространения коронавируса COVID-19 в Москве использовались данные по распространению ви- руса в Китае [11]. Начальное значение нормированного мно- жителя N для Китая подбирали, использовав соотношение (5). Для Китая количество заболевших на конец эпидемии составляло ȳ = 82500. Начальное приближение среднего зна- чения параметра показателя роста популяции брали равным <α> ~ 1,12. Из тренда кривой распределения численности популя- ции инфицированных COVID-19 на протяжении всего пери- ода эпидемии определяли показатели α и значение нор- мировочного множителя N. Интервал времени Δ = tn + 1 - tn , с которым рассчитывали численность популяции инфици- рованных, брали равным 12 часов, а сравнение со стати- стическими данными проводилось через 2Δ, т.е. раз в сут- ки. Значение нормировочного множителя соответствовало N = 760 000, а параметр роста численности популяции инфи- цированных равнялся: α = 1,19 с 3.01.2020 по 30.01.2020; α = 1,104 с 31.01.2020 по 10.02.2020; image N n  1 n     α = 1,18 с 10.02.2020 по 12.02.2020; где yn = y(tn) - численность популяции заболевших в n-й мо- мент времени tn ; λ - коэффициент (показатель) роста числен- ности популяции; N - максимальное значение численности жителей, которые потенциально могут заразиться коронави- русом COVID-19. Этот параметр зависит от целого ряда фак- торов, таких как численность населения, его скученность или плотность, устойчивость к заболеванию, дисциплинирован- ность населения во время карантинных мероприятий и др. Эта модель показывает, что численность популяции за- болевших быстро растет, пока она мала (yn  N) и начинает убывать, когда заболевших становится много. Сделаем заме- ну переменных yn = xnN; α = λ. (2) Тогда уравнение (1) преобразуется к виду xn + 1 = αxn (1 - xn), (3) где переменные xn и параметр α являются безразмерными. При значениях 0 < α ≤ 1 независимо от выбора x1 чис- ленность популяции стремится к нулю [9]. Таким образом, α = 1,119 с 13.02.2020 и до конца эпидемии. На рис. 1 и 2 представлены графики изменения числен- ности популяции инфицированных коронавирусом COVID-19 в Китае и прирост численности в день на период эпидемии. Синие линии - расчет, красные - реальные данные [11]. Из рис. 1 и 2 видно хорошее соответствие расчетных и фактических данных. Это совпадение дало возможность использовать логистические уравнения (1) и (3) для моде- лирования распространения эпидемии коронавируса в ряде европейских и азиатских стран. Для каждой из стран были найдены показатели роста популяции α и потенциальные численности N. В табл. 1 приведены найденные средние показатели ро- ста популяции заболевших до пика эпидемии в ряде евро- пейских и азиатских стран и США. Также был определен нормировочный множитель N для ряда европейских и азиатских стран. Он был найден из со- поставления расчетных и фактических данных. Фактические данные были взяты на сайте официальной статистики [12]. image 100000 Общее число заболевших 80000 60000 40000 20000 Расчет Факт 0 15.01 25.01 04.02 14.02 24.02 05.03 15.03 Дата Рис. 1. Кривая численности инфицированных коронавирусом COVID-19 на протяжении эпидемии 16000 Число новых заболевших 14000 12000 10000 8000 6000 4000 2000 0 image Ряд1 Ряд2 15.01 25.01 04.02 14.02 24.02 05.03 15.03 Дата Рис. 2. Кривая прироста численности инфицированных Таблица 1 Страна Показатель α Страна Показатель α Германия 1,190 Франция 1,149 Португалия 1,180 США 1,148 Испания 1,163 Южная Корея 1,121 Китай 1,160 Швеция 1,088 Италия 1,152 Япония 1,047 Средние показатели роста численности заболевших до пика эпидемии по странам Таблица 2 Соотношение максимально возможного числа потенциально заразившихся вирусом в городе или стране к их общей численности × 104 image 8 7 6 5 4 3 2 1 0 Заболевшие 1 - N = 1 000 000 2 - N = 760 000 3 - N = 500 000 4 - N = 300 000 4 сценария Xend 1 - 75 000 2 - 58 000 3 - 38 750 4 - 25 780 Страна, город N/A ∙ 100% (А - численность города, страны) г. Нью-Йорк, США 11,600 г. Ухань, Китай 4,750 Швеция 5,000 Испания 4,400 Италия 3,660 Португалия 3,000 Южная Корея 1,960 Германия 1,340 США 1,200 Япония 0,550 Китай 0,053 1 2 После проведения анализа соотношения нормировочного 3 множителя к числу жителей Москвы рассмотрены 4 сценария развития распространения эпидемии коронавируса COVID-19. 4 1 сц енари й: N = 1 000 000 чел. (N/A ∙ 100% = 7,1 %); (6.1) 2 сц енари й: N = 760 000 чел. (5,4%); (6.2) 3 сц енари й: N = 500 000 чел. (3,6%); (6.3) 4 сц енари й: N = 300 000 чел. (2,1%). (6.4) 0 10 20 30 40 50 60 70 80 90 Рис. 3. Кривые распределения численности инфицированных по дням эпидемии в соответствии с 4-я сценариями, Москва, РФ На рис. 3 и 4 представлены фактические данные (с 05.03.2020 по 09.04.2020) [13] и расчетные кривые распре- деления численности заболевших и их прироста на протяже- нии всех дней эпидемии (прогноз). 3500 3000 2500 2000 1500 1000 500 image Прирост 1 - N = 1 000 000 2 - N = 760 000 3 - N = 500 000 4 - N = 300 000 24.04 1 22.04 2 19.04 3 19.04 4 Максимум прироста dX X 1 - 3033 39 130 В табл. 1 приведены средние расчетные показатели роста 2 - 2420 30 630 численности заболевших (при которых получено наилучшее 3 - 1626 20 420 согласование с фактическими данными) до пика эпидемии, 4 - 1210 13 636 т.е. те показатели роста, которые определили и время пика, и численность на пике. Из табл. 1 следует, что наиболее низ- кие показатели роста численности заболевших у стран Шве- ция и Япония. Наверное, благодаря этим показателям они и не вводили жесткие карантинные меры. На развитие сцена- рия эпидемии влияет нормировочный множитель, характе- ризующий максимальное значение числа жителей, которые потенциально могут быть инфицированы. Это не означает, что все N жителей города или страны заболеют, заболеет столько, сколько соответствует показателю роста α. Нормиро- моделирования прогноза распространения эпидемии в ев- ропейских странах: Испания, Италия, Франция, Португалия, Швеция; азиатских странах: Южная Корея, Япония и США. 0 0 20 40 60 80 100 Рис. 4. Кривые прироста числа заболевших по дням эпидемии в соответствии с 4-я сценариями, Москва, РФ Значения параметров показателей роста численности по- пуляции заболевших представлены ниже. Сц енарий 1: N = 1 000 000 чел. α = 1,14 с 5 по 21 марта; α = 1,111 с 22 марта по 1 апреля; α = 1,081 с 1 апреля. Сц енарий 2: N = 760 000 чел. α = 1,14 с 5 по 21 марта; α = 1,111 с 22 марта по 1 апреля; α = 1,083 с 1 апреля. Сц енарий 3: N = 500 000 чел. α = 1,14 с 5 по 21 марта; α = 1,113 с 22 марта по 1 апреля; α = 1,084 с 1 апреля. Сц енарий 4: N = 300 000 чел. α = 1,14 с 5 по 21 марта; α = 1,108 с 22 марта по 1 апреля; α = 1,094 с 1 апреля. ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ Наличие хорошего соответствия между фактическими данными распространения эпидемии коронавируса COVID-19 в Китае и расчетными результатами на основании логистиче- ского уравнения (1) дало возможность его использовать для вочный множитель N зависит от иммунитета нации к вирусу, условий проживания жителей (скученность и т.п.), менталите- та и т.д. и рассчитывается на основании согласования расчет- ных и фактических данных. В табл. 2 представлено отноше- ние N к численности населения A выбранных городов и стран (N/A ∙ 100%). Для определения нормировочного множителя N для Москвы была проанализирована табл. 2 и предложе- но четыре возможных сценария распространения COVID-19 в Москве (6.1)-(6.4) с соотношением (N/A) от 2 до 7%. Первый сценарий (6.1) - самый тяжелый, назовем его Нью- Йоркским, второй - Уханьским (нормировочный множитель совпадает с нормировочным множителем, взятым для Уханя). Последний четвертый сценарий - самый легкий, назовем Из- раильским (нормировочный множитель как у Израиля). В табл. 3 для 4-х сценариев показано время пика эпиде- мии коронавируса в Москве, численность на пике, прирост числа заболевших на пике и численность заболевших к кон- цу эпидемии. В последнем столбце указано расчетное время окончания эпидемии (когда заболевает только несколько че- ловек в день). Показатель роста численности популяции заболевших α определяется на ранних стадиях распространения эпи- демии, когда в уравнении (1) преобладает линейный член, а нелинейным можно пренебречь: image xn  1  . (7) xn Тогда заболевание нарастает по экспоненте, которая в логарифмической шкале представляет собой прямую. На- клон прямой задается показателем роста. Если показатель роста изменяется, то изменяется и наклон прямой. На рис. 5 в логарифмической шкале представлен рост числа заболев- ших по дням (синяя кривая - статистические данные, крас- ная - расчеты по модели для сценария с N = 300 000). Таблица 3 Параметры пиков эпидемии для 4-х сценариев, Москва, РФ № сценария N Время пика Численность на пике, чел. Прирост числа заболевших на пике, чел. Число заболевших в конце, чел. Время окончания эпидемии 1 1 ∙ 106 24,04 39 130 3033 75 000 Середина июня 2 7.6 ∙ 105 22,04 30 630 2420 58 000 Начало июня 3 5 ∙ 105 19,04 20 420 1626 38 750 Начало июня 4 3 ∙ 105 15,04 13 636 1210 25 780 Конец мая Мы видим, что весь отрезок наблюдения можно прибли- женно разбить на три интервала с разными показателями α, в которых тренд роста заболевших имеет разный наклон. Значения показателя α представлены в табл. 4. Когда накапливается большое количество заболевших, в уравнение (1) начинает играть роль нелинейный член, и от- ношение имеет вид x image n  1  1 - x . (8) x 105 image n n Заболевшие N = 300 000 104 103 01.04 09.04 α = 1,094 с 01.04 Все рассмотренные сценарии для Москвы хорошо описывают статистические данные на сегодняшний день. По ка- кому из сценариев пойдет развитие эпидемии станет ясно ближе к пику. Следует отметить, что прогноз осуществлен со значени- 102 101 100 22.03 α = 1,143 с 5.03 по 21.03 α = 1,108 с 21.03 по 01.04 Данные ем показателя не меняющимся с 1 апреля. А значение этого параметра коррелируется с количеством контактирующих людей. Мы видим, что принятые меры по самоизоляции на- селения произвели эффект, и во всех сценариях показатель роста уменьшился. ЗАКЛЮЧЕНИЕ 0 5 10 15 20 25 30 35 40 Рис. 5. Определение показателя численности популяции инфицированных коронавирусом COVID-19, Москва, РФ Таблица 4 Значение показателя роста популяции Значение показателя α Время его удержания 1,143 05.03-21.03 1,108 22.03-01.04 1,094 02.04 - настоящее время Из табл. 4 видно, что значения показателя численности заболевших падают с введением мер по изоляции и умень- шению контактов. Показано использование дискретного логистического уравнения для моделирования распространения корона- вируса COVID-19 в Москве. Определены показатели роста численности популяции инфицированных коронавирусом COVID-19 для 4-х сценариев развития эпидемии в Москве. Для «легкого» варианта - рассчитанное время пика прихо- дится на 15.04 с приростом числа заболевших на пике 1210 чел. и общей численностью 13 636 чел. Для «тяжелого» вари- анта время пика - 24.04 с приростом числа заболевших 3033 чел. и общей численностью на пике 38 130 чел. Осуществлен прогноз численности заболевших для сценариев: «легко- го» - 25 780 чел., «тяжелого» - 75 000 чел. Какое значение будет принимать нормировочный множитель N (характери- зующий максимальную численность потенциально возмож- ных инфицированных вирусом COVID-19 жителей Москвы) станет ясно ближе к времени прохождения пика эпидемии, когда заболевших будет достаточно много и начнет влиять нелинейность логистического уравнения.

About the authors

Eleonora M. Koltsova

Mendeleev University of Chemical Technology of Russia

Email: kolts@muctr.ru
Doctor of Engineering, Professor; Head of Department IСT Moscow, Russian Federation

Elena S. Kurkina

Mendeleev University of Chemical Technology of Russia; Lomonosow Moscow State University

Email: e.kurkina@rambler.ru
Doctor of Physics and Mathematics, Associate Professor; professor of Department IСT; leading researcher of Department BMK Moscow, Russian Federation

Aleksey M. Vasetsky

Mendeleev University of Chemical Technology of Russia

Email: amvas@muctr.ru
senior lecturer of Department IСT Moscow, Russian Federation

References

Supplementary files