Журнал вычислительной математики и математической физики, 2021, T. 61, № 3, стр. 400-412

1. ВВЕДЕНИЕ

Математические модели распространения инфекции имеют особое значение в случае глобальных пандемий. Несмотря на то что они не могут воспроизвести полную картину пандемии, они обеспечивают качественные прогнозы ее развития и выявляют конкретные динамические закономерности, которые могут быть полезны при оптимальном планировании социальных профилактических мер, которые позволят в определенной степени контролировать ситуацию. Эти меры являются предметом жарких споров, поскольку они существенно меняют повседневную жизнь людей и имеют большие экономические последствия. Отметим, например, что эффективность социального дистанцирования для снижения пика заболеваемости была продемонстрирована еще во время пандемии испанского гриппа в 1918–1919 гг. [1]. Новые инструменты прогнозирования, основанные на математическом моделировании, обеспечивают несомненное преимущество в облегчении процесса адаптации и корректировки, способствуя эффективному управлению ресурсами на индивидуальном и институциональном уровнях.

Попытки математического моделирования инфекций имеют более чем столетнюю историю [2]–[4]. Позже Кермак и МакКендрик [5] рассмотрели развитие эпидемии в закрытой однородной популяции, предполагая, что полный иммунитет обеспечивается однократным заражением, и что в сам момент заражения человек не заразен. С этими предположениями проблема была в конечном итоге сведена к интегродифференциальному уравнению типа Вольтерра, анализ которого позволил сделать такие выводы, как существование пороговой плотности населения, ниже которой не может возникнуть эпидемия, решающая роль небольших изменений скорости заражения, и конец эпидемии до того, как уязвимое население будет исчерпано. Модель была обобщена [6] с учетом эффекта от постоянного притока новых восприимчивых особей в популяцию (рождение, иммиграция и т.д.). Аналогичные результаты были получены при передаче инфекции через промежуточного хозяина [5].

Естественно, пандемия COVID-19 вызвала всплеск интереса к разработке инструментов прогнозирования, основанных на математическом моделировании. Конкретные цели и подходы разных исследовательских групп существенно различаются. Например, Uhlig и соавт. [7], основываясь на детерминированных компартментных моделях, сочетающих эпидемиологический, статистический и нейросетевой подходы, предложили метод эмпирического нисходящего (top-down) моделирования для обеспечения прогнозов эпидемий и расчетов рисков для (локальных) вспышек. Основываясь на первоначальных результатах, авторы [7] предполагают, что статистические системы предлагаемого типа будут использованы для управления автоматическими веб-платформами с целью демократизации распространения результатов прогнозов.

Для оперативных оценок желательны более простые режимы. Самый простой подход – это подгонка графиков к имеющимся данным. В [8] предполагалось, что количество $P(t)$ совокупных диагностированных положительных случаев COVID-19 описывается функцией ошибок. Это верно, если ежедневный прирост новых случаев $P{\kern 1pt} '(t)$ может быть описан гауссовым распределением, что обычно не так [9]. Köhler-Rieper и соавт. [9] предложили детерминированную модель прогноза для эпидемии COVID-19, в которой динамика модели выражается одной прогностической переменной $\kappa (t)$, входящей в интегродифференциальное уравнение в виде коэффициента, зависящего от времени. Модель [9] имеет сходство с классическими компартментными моделями, такими как SIR [10], а переменная $\kappa (t)$ может быть интерпретирована как эффективное число воспроизводства случаев инфицирования. Авторы [9] считают это принципиальным преимуществом, поскольку модель сформулирована с использованием наиболее достоверных статистических данных, а именно, количества кумулятивных диагностированных положительных случаев COVID-19. Они применили модель к более чем 15 странам, и результаты доступны через веб-платформу [11].

В качестве альтернативы интегродифференциальному уравнению в данной статье мы используем подход, основанный на дифференциальных уравнениях с запаздыванием [12], частном случае функционально-дифференциальных уравнений. Такие модели задержки следуют стратегии моделирования процесса удаления не с помощью отдельной переменной, а со сдвигом по времени в функции, описывающей количество кумулятивных случаев. Dell’Anna [13] представил пример применения модели задержки по времени к проблеме COVID-19. Можно показать эквивалентность интегродифференциальных и дифференциальных уравнений с запаздыванием. Например, авторы [9] показывают, что основное уравнение их модели идентично дифференциальному уравнению с запаздыванием в [7] при соответствующем выборе весовой функции интегрирования.

Наше модельное уравнение аналогично общей модели из [5], [13], однако последующий анализ в [13] проводится с использованием упрощенного уравнения без учета так называемого логистического фактора (см. ниже), т.е. в частном случае пренебрежимо малого количества инфицированных лиц по сравнению с численностью популяции, что реально в глобальном смысле, но сомнительно для небольших групп населения

2. ОСНОВНЫЕ ОПРЕДЕЛЕНИЯ

${{N}_{{max}}}$ – количество членов популяции.

$N(t)$ – количество инфицированных в момент времени $t$.

$\tau $ – время возможного распространения инфекции одним вирусоносителем. Это может быть временем естественной длительности болезни или временем от заболевания до изоляции вирусоносителя от общества. В настоящем рассмотрении $\tau $ – параметр модели.

$N(t - \tau )$ – количество больных, но уже не вирусоносителей. Разумеется, при $t - \tau < 0$ величина $N(t - \tau ) = 0$.

$N(t) - N(t - \tau )$ – количество вирусоносителей в момент времени $t$.

${{N}_{{max}}} - N(t)$ – количество не зараженных членов популяции.

$(N(t) - N(t - \tau )){\text{/}}{{N}_{{max}}}$ – плотность вирусоносителей в обществе.

$\alpha $ – вероятность инфицирования здорового члена популяции при контакте с инфицированным в единицу времени, например за день. Эта вероятность определяется как вероятность инфицирования при одном контакте с носителем инфекции, умноженная на количество контактов человека со всеми членами популяции в единицу времени. В настоящем рассмотрении $\alpha $ – параметр модели.

$P = \alpha \Delta t(N(t) - N(t - \tau )){\text{/}}{{N}_{{max}}}$ – вероятность инфицирования за единицу времени $\Delta t$, например за 1 день, для одного члена популяции при заданной плотности носителей инфекции.

Тогда количество заболевших за единицу времени

Видно, что в этой модели можно использовать не количество членов популяции $N(t)$, а понятие плотности, $N(t){\text{/}}{{N}_{{max}}} = x(t)$, что делает ее пригодной для описания любого сообщества (страна, город и т.п.). В этом случае уравнение (3) для плотности $x(t)$ принимает вид

В уравнении (4) есть два параметра модели $\alpha $ и $\tau $. Простым введением “универсальной шкалы времени” $\zeta $ в единицах $\tau $: $t = \zeta \tau $, уравнение (4) может быть переписано в виде

Сделаем еще одно замечание. Уравнение (5), равно как и уравнения (4) и (3), в нашей постановке имеет аналитическое решение. Чтобы увидеть это, введем для наглядности обозначения:

Аналитические решения, приведенные выше, пригодны для полного анализа развития пандемии в рамках описанной модели, но для этого необходима работа со сложными и длинными в записи выражениями. Единственная польза от этого рассмотрения, без анализа общего решения, состоит в том, что при непрерывности функции $x(t)$, задаваемой разбиением на решениях в интервалах, длительностью $\tau $, производная $dx(t){\text{/}}dt$ тоже будет непрерывной функцией за исключением одной точки $t = \tau $. В этой точке $x(t - \tau ) = x(0)$ не может быть равной $0$ (при таких начальных условиях существует только тривиальное решение $x = 0$) и производная, определяемая уравнениями (4) или (5), всегда будет иметь разрыв.

Обговорив возможности аналитического решения, далее мы, на этом этапе рассмотрения, использовать его не будем из-за его сложности. Будем исследовать достаточно очевидные свойства решений в линейном приближении или проводить анализ численных решений уравнений (4) и (5).

3. ПРЕДВАРИТЕЛЬНОЕ ИССЛЕДОВАНИЕ

3.1. Начальная стадия

Видно, что уравнение (4) нелинейно и нелокально из-за слагаемого, сдвинутого по времени. Тем не менее на временах, когда заражение не охватило все общество и $x(t) \ll 1$, уравнение (4) может быть исследовано в линейном приближении:

(6)

Начнем исследование с начала инфицирования, когда $x(t - \tau ) = 0$. В этом случае уравнение (6) имеет экспоненциальное решение:

(7)

$x(t) = x(0){{e}^{{\alpha t}}}$

с количеством относительного прироста инфицированных (например, за 1 день, $dt \approx \Delta t = 1$ день)

(8)

$\frac{{dx}}{x} = \alpha dt\quad {\text{или}}\quad \frac{{\Delta x(t)}}{{x(t)}} = \alpha ,$

т.е. коэффициент $\alpha $ определяет относительный рост инфицированных в начале пандемии. Его значения достигают $0.3{\kern 1pt} - {\kern 1pt} 0.2$ (см. любую страну по данным института Хопкинса). С вводом режима (само)изоляции значения коэффициента $\alpha $ понижаются и становятся близкими к $0.15{\kern 1pt} - {\kern 1pt} 0.1$ (опять же – данные института Хопкинса). Для понимания дальнейшего изменения темпов роста инфицированных необходимо уже анализировать уравнение (6) вместе с запаздывающим слагаемым. Обратим внимание, что решение

(9)

$x = {{x}_{0}}{{e}^{{\beta t}}}$

удовлетворяет уравнению (6), если $\beta $ удовлетворяет трансцендентному уравнению:

(10)

$\beta = \alpha (1 - {{e}^{{ - \beta \tau }}}),$

которое легко переписывается в уравнение для $z = \beta \tau $:

(11)

$z = \alpha \tau (1 - {{e}^{{ - z}}}).$

Поэтому, вообще говоря, относительное количество инфицированных за 1 день через несколько недель после начала пандемии будет определяться параметром $\beta $:

(12)

$\frac{{\Delta x(t)}}{{x(t)}} = \beta ,$

который, в свою очередь, определяется уравнением (11). При $\alpha \tau \gg 1$ уравнения (10) и (11) имеют решения $\beta \approx \alpha $, т.е. решения уравнений (7) и (9) совпадают. Заметим, что при $\alpha \tau = 3$, $\beta \approx 0.94\alpha $, а при $\alpha \tau = 2.6$, $\beta \approx 0.9\alpha $, т.е. $\beta $ практически равно $\alpha $ (в рамках точности данных, приведенных в Википедии) уже при значении параметра $\alpha \tau = 2.5$. Второе очевидное решение уравнения (10) $\beta = 0$ порождает решение в виде константы и ее вкладом на экспоненциально растущем участке слагаемым можно пренебречь. За исключением случая, когда у нас имеется только одно решение уравнений (10) и (11). Это происходит, когда $\alpha \tau \leqslant 1$. В этом случае решение $\beta = 0$ единственное и возникает возможность прекращения пандемии на стадии вне полного инфицирования.

4. ПРИМЕРЫ ЧИСЛЕННЫХ РЕШЕНИЙ

4.1. Решения при $\alpha \tau > 1$

Для численного исследования уравнение (4) решалось методом Эйлера с шагом в 1 день. На фиг. 1 приведены численные решения для ситуации “свободный полет”, когда $\alpha \tau = 2.8 > 1$. Может иллюстрировать политику полной иммунизации общества вследствие его тотального инфицирования. Начальные условия выбраны равными одной миллионной от населения (320 для США, 140 для России, 12 для Москвы, 2.5 для Алматы и т.д.). Время инфицирования выбрано равным 14 дням. Скорее всего, нужно брать 20–30 дней.

Фиг. 1.

(а) Изменение плотности инфицированных $x(t)$ со временем $t$ в днях, (б) абсолютный прирост плотности инфицированных $dx{\text{/}}dt$, (в) относительный прирост плотности инфицированных $(dx{\text{/}}dt){\text{/}}x$ при $\alpha = 0.2$, $\tau = 14$ дней, $\alpha \tau = 2.8$, $x(0) = {{10}^{{ - 6}}}$.

Результат: длительность пандемии около 100 дней. Пик заболевания приходится на 80–90-е дни и достигает значений 40 человек на 1000 населения. Поскольку такого количества больничных коек нигде нет (согласно сайту ВОЗ [14] количество больничных коек не превышает 8 на 1000 населения даже для развитых стран), то эта схема будет сопровождаться большим количеством отказов в лечении (т.е. смертей). Интересно, что относительный прирост сопровождается скачком при переходе на 14-й день решения (7) в решение (9), т.е. небольшим отличием $\alpha $ от $\beta $ (см. выше). О разрыве производных в этой точке говорилось выше.

Заметим, что при $\alpha \tau = 2.8$ уравнение (14) дает, что пик заболевания будет при достижении ${{x}_{m}} = 0.64$, т.е., согласно фиг. 1а, на 86-й день развития эпидемии. Этот вывод подтверждается фиг. 1б, где пик заболевания приходится на 82-й день.

4.2. Решения при $\alpha \tau < 1$

Для демонстрации режима угасания пандемии при $\alpha \tau < 1$, когда существует единственное решение уравнения (10): $\beta = 0$, ниже приводится фиг. 2 с расчетами уравнения (4) при $\alpha = 0.05$, $\tau = 14$ дней, т.е. $\alpha \tau = 0.7 < 1$. Начальное значение было выбрано $x(0) = {{10}^{{ - 5}}}$, т.е. в 10 раз больше, чем для расчетов, приведенных на фиг. 1. Это было сделано для удобства анализа, поскольку рост инфицированных небольшой. Видно, что рост инфицированных (примерно в 3.5 раза) выходит на константу, не достигая единицы, т.е. возможен режим остановки пандемии без полного инфицирования популяции. Для этого произведение $\alpha \tau $ должно быть меньше 1, как и обсуждалось выше. В этом режиме производная имеет заметный скачок, поскольку $\alpha $ и $\beta $ значительно различаются. Для этого примера продолжительность пандемии около 100 дней и близка ко времени развития пандемии, приведенной на фиг. 1. Это случайность. Можно выбрать вероятность инфицирования больше, например, $\alpha = 0.07$ и получить значительно большее время развития пандемии (около 1500 дней) в режиме $\alpha \tau < 1$. Количество инфицированных при этом увеличится примерно в 50 раз, т.е. решение, по-прежнему, будет выходить на константу, значительно отличающуюся от 1, но за очень большие времена.

Фиг. 2.

(а) Изменение плотности инфицированных $x(t)$ со временем $t$ в днях, (б) абсолютный прирост плотности инфицированных $dx{\text{/}}dt$, (в) относительный прирост плотности инфицированных $(dx{\text{/}}dt){\text{/}}x$ при $\alpha = 0.05$, $\tau = 14$ дней, $\alpha \tau = 0.7$, $x(0) = {{10}^{{ - 5}}}$.

4.3. Зависимость решений от значения параметра $\alpha \tau $

Рассмотренные примеры численных решений не противоречат анализу уравнений (4) в линейном приближении. Более того, подтверждают основные выводы. Тем не менее остается неясным выход на константу при больших временах для $\alpha \tau > 1$. В этом случае для анализа нельзя использовать линейное приближение и нужно проводить анализ численно. Поскольку уравнение (5) определяется только одной константой $\alpha \tau $, то ниже приводятся численные решения при разных значениях этой константы. Правильная схема отображения плотности инфицированных в зависимости от значения $\zeta $ не является наглядной. Поэтому ниже приведена ось времени, совпадающая с днями при $\tau = 14$ дней. При других значениях $\tau $ множитель $\tau {\text{/}}14$ является масштабирующим. Например, 100 дней при $\tau = 14$ перейдут в 200 дней при $\tau = 28$. Во всех расчетах $x(0) = {{10}^{{ - 6}}}$.

Фиг. 1 при $\alpha \tau = 0.5,\;0.9,\;0.99$ показывают динамику плотности инфицирования и ее производную по времени $t$ при $\alpha \tau < 1$. В этом случае пандемия прекращается при малом количестве инфицированных (максимум 0.01% популяции), но время ее развития достаточно велико и увеличивается при приближении $\alpha \tau $ к 1. На этих фигурах при $\alpha \tau $ = 1.5, 2, 3, 5 для сравнения показаны кинетики нарастания инфицированных, для которых доля инфицированных меняется от 60 до 100% за достаточно короткие времена развития инфекции. На фиг. 3 также приведены решения при $\alpha \tau = 1.01,\;1.1$ в случае небольшого отличия $\alpha \tau $ от 1. Так, при $\alpha \tau = 1.01$ решение выходит на константу (2% инфицированного населения) в течение длительного периода (годы). При $\alpha \tau = 1.1$ срок выхода на плато короче, но количество инфицированных составит 18%. Из фиг. 3б видно, что при $\alpha \tau > 1$ имеет место возрастание производной $dx{\text{/}}dt$ на некотором интервале при $t > \tau $, чего не наблюдается при $\alpha \tau \leqslant 1$. Заметим, что в настоящий момент. (Статья подготовлена в мае 2020 г.) Текущая статистика COVID-19 за более поздние периоды (см., например, [15], [16]), качественно согласуется с нашими прогнозами) количество инфицированных в России составляет 0.13%, Германии – 0.20% в США – 0.36%, т.е. в этих странах точно реализуется сценарий $\alpha \tau > 1$.

Фиг. 3.

(а) Численные решения $x = x(t)$ уравнения (5) и (б) их производные $dx{\text{/}}dt$ в зависимости от $t$ в единицах $\tau /14$ при разных значениях $\alpha \tau $: 0.5, 0.9, 0.99, 1.01, 1.1, 1.5, 2, 3, 5.

4.4. Волны эпидемии

Приведенные выше рассуждения показывают, что при заданной константе $\alpha \tau $ пик заболевания определяется уравнением (13) и является единственным. Однако при возможности управления параметром $\alpha \tau $ можно добиться двух и даже трех пиков. Для демонстрации этого утверждения введем управление параметром $\alpha $ с резким изменением его значения:

(18)

$\alpha = {{\alpha }_{0}}(1 + \varepsilon - \varepsilon {{e}^{{ - {{{(t/T)}}^{6}}}}}).$

Для выбора значения $\varepsilon $ можно опираться на значения эффективной константы $\alpha \tau (1 - {{x}_{0}})$, которая должна быть больше 1 и асимптотического значения ${{x}_{\infty }}$, например, из фиг. 3 при ${{x}_{0}} = {{x}_{\infty }}$. Ниже приводятся (см. фиг. 4), иллюстрирующие расчеты с появлением двух и даже трех пиков. Видно, что при выбранном значении ${{\alpha }_{0}}\tau $ вблизи 1 небольшой скачок $\alpha $ на 20% (фиг. 4а и 4б) приводит к появлению второго пика (фиг. 4б), а для того чтобы они были наблюдаемы, время переключения параметра $\alpha $ составляет $T = 1500$ дней (сравни фиг. 4а и 4б). Для появления 3-го пика скачок $\alpha $ должен быть больше. В приведенном примере (фиг. 4в) $\alpha $ увеличивается в 3 раза.

Фиг. 4.

Абсолютный прирост $dx{\text{/}}dt$ плотности инфицированных в зависимоcти от $t$ в единицах $\tau $: (а) при $T = 1100$ дней и $\alpha = {{\alpha }_{0}}(1.2 - 0.2exp( - {{(t{\text{/}}T)}^{6}}))$, (б) при $T = 1500$ дней и изменении $\alpha = {{\alpha }_{0}}(1.2 - 0.2exp( - {{(t{\text{/}}T)}^{6}}))$, (в) при $T = 1500$ дней и $\alpha = {{\alpha }_{0}}(3 - 2exp( - {{(t{\text{/}}T)}^{6}}))$ при ${{\alpha }_{0}}\tau = 1.1$, $\tau = 14$ дней, $x(0) = {{10}^{{ - 6}}}$.

Выбор переходной функции в форме (18) является достаточно произвольным, и можно исследовать функции с более быстрой степенью переключения значения $\alpha \tau $, сжимая, таким образом, шкалу времени для возникновения новых пиков заболевания. Именно за это поведение отвечает реакция общества на пандемию.

Таким образом, ответ на вопрос: “Можно ли ожидать двух или более пиков скоростей заболевания?”, можно сформулировать в форме: “Количество пиков может достигать, по крайней мере, трех, и определяется амплитудой и скоростью изменения параметра $\alpha \tau $.”

4.5. Транспорт инфекции

Рассматривая развитие пандемии в рамках уравнений (4), (5), т.е. в схеме “инфицированный больной заражает окружающих”, мы не учитывали приток инфицированных больных извне. И если для больших сообществ этот механизм роста больных может оказаться малозначимым, то для небольших городов приток больных извне может оказывать существенное влияние на динамику развития инфекции. Для того, чтобы включить поток инфицированных в уравнения (4) или (5), можно просто добавить слагаемое $j$ к производной $x(t)$:

(19)

$\frac{{dx}}{{dt}} = (1 - x(t))(x(t) - x(t - \tau )) + j,$

где $j$ имеет смысл относительной доли инфицированных больных, прибывающих в сообщество в единицу времени, например, за день. В общем случае поток инфицированных больных $j$ является функцией времени $t$ и может быть легко изменен усилиями общества, вплоть до полной изоляции ($j = 0$). При рассмотрении схемы развития инфекции по уравнению (19) начальное условие $x(t = 0)$ может быть нулевым, убирая описанный выше скачок производной при $t = \tau $. Используемые выше схемы анализа решений уравнения (4) пригодны и для анализа уравнения (19). Мы не будем их здесь повторять, а приведем пример численного решения уравнения (19), отображенного на фиг. 5, при относительно небольших значениях $\alpha \tau $ и $j$. Выбранное значение $j = {{10}^{{ - 6}}}$ означает, что в город с населением сто тысяч человек въезжает один инфицированный человек за десять дней. Как видно из фиг. 5, даже такой незначительный приток инфицированных людей приводит к вспышке заболевания, протекающей быстрее, чем при задании начального условия на количество больных. Это можно увидеть, сравнивая фиг. 5 с кривой при $\alpha \tau = 1.1$ на фиг. 3а.

Фиг. 5.

(а) Плотность $x(t)$ инфицированных в зависимости от времени $t$ в единицах $\tau $, и (б) абсолютный прирост $dx{\text{/}}dt$ плотности инфицированных при $\alpha \tau = 1.1$, $\tau = 14$ дней, $x(0) = 0$, $j = {{10}^{{ - 6}}}$.

Характерной особенностью динамики заболеваний при внешнем притоке носителей инфекции является “затянутый хвост”, который, в принципе, позволяет увидеть источник заболеваний во внешнем притоке носителей инфекции.

5. СРАВНЕНИЕ НАСТОЯЩЕЙ МОДЕЛИ С МОДЕЛЬЮ SIR

Как уже отмечалось выше, пандемия “испанки” спровоцировала необходимость создания прогноза пандемий на основе хоть и упрощенных, но математических моделей. В частности, именно тогда была создана до сих пор используемая модель SIR [5]. В настоящее время эта модель значительно расширена, но ее уравнения до сих пор остаются базовыми, именно с ними мы будем сравнивать настоящую модель. Название модели является аббревиатурой трех переменных модели: $S$ – Susceptible (количество людей, способных заразиться), $I$ – Infectious (количество носителей инфекции), $R$ – removed (количество людей, не способных заразиться (иммунные или умершие)). Система дифференциальных уравнений, восходящая еще к [5], записывается в виде

Таким образом, уравнение нашей модели не противоречит модели SIR, являясь ее реализацией при функциональной связи между заболевшими и носителями инфекции $I = N(t) - N(t - \tau )$.

6. АНАЛИЗ ТЕКУЩЕЙ СИТУАЦИИ

В настоящий момент считается, что антитела на инфекцию возникают практически у всех заболевших на 20–24-й день после инфицирования. Поэтому в этом анализе будем полагать, что время инфицирования $\tau $ составляет 21 день в социуме, где противоинфекционные меры общества не позволяют его снизить радикально (“свободный полет”). С учетом оценочной вероятности инфицирования в режиме (само)изоляции $\alpha = 0.1{\kern 1pt} - {\kern 1pt} 0.15$ получаем значения $\alpha \tau = 2.2{\kern 1pt} - {\kern 1pt} 3.15$, которые дают (см. фиг. 3 при $\alpha \tau = 2,\;3$) время инфицирования популяции примерно 150 дней с естественной датой завершения пандемии (февраль, март, апрель, май, июнь) конец июня–июль практически для всех стран после вспышки в Китае. Если же общество сократит время инфицирования $\tau $ в 2–3 раза, то длительность пандемии увеличится и будет определяться единицей времени “год”, с бонусом в меньшее количество заболевших.

Разумеется, действия общества по уменьшению $\alpha \tau $ приведут к необходимости решения уравнений (4) при переменном значении $\alpha \tau $, и приведенные численные решения не будут отображать ситуацию в полной мере. В этом случае с целью прогноза можно решать уравнения “в текущий момент времени”. Для такой задачи необходимо будет знание $\alpha \tau $ в текущий момент времени. Заметим, что действия общества приблизили значения $\alpha \tau $ к 1 и наблюдаемое увеличение инфицированных определяется параметром $\beta $, т.е. решением уравнения (10). В случае $\alpha \tau $, близких к 1, можно указать линейную связь между этими параметрами:

На фиг. 6 представлены данные из Wikipedia по числу инфицированных COVID-19 в Москве, и вычисленный на основе этих данных согласно уравнению (2) параметр модели $\alpha \tau \equiv \alpha (t)\tau $:

Фиг. 6.

(а) Число инфицированных COVID-19 $N(t)$, (б) суточный прирост $\Delta N(t)$ инфицированных cо 2 марта по 13 августа в Москве по доступным данным из Википедии) и (в) зависимость параметра $\alpha \tau \equiv \alpha (t)\tau $ c 1 апреля по 13 августа при $\tau = 21$ дней.

Относительно равномерное распределение зараженных по Москве к началу эпидемии и медленное распространение, при котором человек, зараженный в одном конце Москвы, может заразить в другом ее конце, а также условия (законы), не зависящие от района Москвы, позволяют пренебречь очаговостью и рассматривать весь мегаполис как одну популяцию. Большое число заболевших минимизируют вероятностную составляющую ошибки, что видно на фиг. 6, 7, где амплитуда быстрых колебаний дневного прироста заболевших и произведения $\alpha \tau $ со временем уменьшается.

Фиг. 7.

(а) Зависимость параметра $\alpha \tau \equiv \alpha (t)\tau $ c 1 апреля по 13 августа для Москвы при $\tau = 7,\;14,\;21,\;28$ дней, (б) число инфицированных $N(t)$ и (в) суточный прирост инфицированных $\Delta N(t)$, полученные с помощью экстраполяции на 2021–2024 гг.

Для иллюстрации на фиг. 7 представлены параметр модели $\alpha \tau \equiv \alpha (t)\tau $ при $\tau = 7,\;14,\;21,\;28$ и прогноз числа $N(t)$ и суточного прироста инфицированных $\Delta N(t)$, полученных, используя среднее значение $\alpha \tau \equiv \alpha (t)\tau $ с 31 июля по 13 августа. При всех $\alpha \tau $ наблюдается возрастание суточного прироста инфицированных $\Delta N(t)$. При этом при $\tau < 14$ cкачок $\alpha \tau $ в окрестности 1 мая становится необъяснимым.

7. ЗАКЛЮЧЕНИЕ. НЕДОСТАТКИ МОДЕЛИ

Предложенная двухпараметрическая модель развития инфекции в виде обыкновенного нелинейного дифференциального уравнения первого порядка с запаздыванием по времени по сути является редуцированной моделью SIR при функциональной связи между заболевшими и носителями инфекции $I = N(t) - N(t - \tau )$. Такая редукция сохраняет оптимальный баланс между адекватностью описания пандемии в модели SIR и простотой практических оценок. При этом предложенная модель позволяет решать как прямую задачу: при заданных параметрах $\tau $ и $\alpha $ находить зависимость плотности инфицированных $x(t)$ от времени $t$, так и обратную задачу: при заданной зависимости плотности инфицированных от времени определить зависимость параметров модели от времени. Это дает возможность быстрого прогноза развития инфекции по предшествующей информации о статистике заболевания (см. фиг. 7).

Отметим, что при заданной константе $\alpha \tau $ пик заболевания определяется уравнением (13), и является единственным. Однако при возможности резкого изменения параметра $\alpha \tau $ можно добиться двух и даже трех пиков, что продемонстрировано на фиг. 4. Такое изменение $\alpha \tau $ за счет изменения $\alpha $ может носить сезонную зависимость восприимчивости к болезни от времени года или зависимость от количества контактов в рабочие и праздничные дни, и дни каникул. В частности, для России минимальное значение $\alpha \tau $ приходилось на весенние и летние каникулы и можно будет видеть в осенние и зимние каникулы. Изменения величины $\alpha \tau $ за счет параметра $\tau $ в большей степени определяются реакцией общества на заболевание. Таким образом, ответ на вопрос: “Можно ли ожидать двух или более пиков скоростей заболевания?”, можно сформулировать в форме: “Количество пиков может достигать, по крайней мере, трех, и определяется амплитудой и скоростью изменения параметра $\alpha \tau $.” Текущая информация, подтверждающая многопиковое распространение COVID-19 в мегаполисах, представлена, например, в сообщениях РИА НОВОСТИ [17].

Рассматривая развитие пандемии в рамках уравнений (4), (5), т.е. в схеме “инфицированный больной заражает окружающих”, мы не учитывали приток инфицированных больных извне. И если для больших сообществ этот механизм роста больных может оказаться малозначимым, то для небольших городов приток больных извне может оказывать существенное влияние на динамику развития инфекции. Для того чтобы включить поток инфицированных в уравнения (4) или (5), достаточно добавить слагаемое $j$ – относительную долю инфицированных больных, прибывающих в сообщество в единицу времени, к производной $x(t)$ и получить в результате неоднородное уравнение (19). В общем случае поток инфицированных больных $j$ является функцией времени $t$ и может быть легко изменен усилиями общества, вплоть до полной изоляции ($j = 0$). Фиг. 5 показывает, что даже незначительный приток инфицированных людей приводит к вспышке заболевания, протекающей быстрее, чем при задании начального условия на количество больных.

К недостаткам модели можно отнести описание в среднем. Примеры проявления некорректности такого описания могут проявляться в небольших социумах, где подпитка инфекции со стороны может оказать большее значение, чем механизмы развития, представленные уравнениями (4). Кроме того, описание не зависит от “координат”, т.е. не отражает пространственного распределения инфекции. Поэтому вне модели остались последующие вспышки из-за “выравнивания плотности инфицированных” по регионам. Остальные недостатки, связанные с количеством параметров модели (например, распределение населения на разные группы со своими значениями $\alpha \tau $), скорее всего не являются значимыми при описании популяции “в целом”. Существенная особенность COVID-19, которая делает его в некоторой степени уникальным среди других инфекций, заключается в том, что тяжесть заболевания сильно различается у отдельных пациентов от легкого до смертельного исхода. Более того, в процессе развития болезни есть особая (бифуркационная) точка, в которой у некоторых пациентов неожиданно возникают фатальные легочные осложнения. Эти особенности имеют большое практическое значение, поскольку в простых случаях статистические данные сильно недооцениваются (пациенты не обращаются за медицинской помощью и не принимаются во внимание), а в потенциально опасных случаях требуется особая немедленная помощь. В нынешнюю упрощенную модель эти особенности, природа которых еще далеко не ясна [18], не были включены.

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

Авторы благодарны А. Кукетаеву за идею взаимосвязи данной модели и модели SIR.