Основное содержание: модель описывает переходы между численностью индивидов, разделенных на три группы: восприимчивые (Susceptible), которые при контакте с инфицированными (Infected) с вероятностью $k$ заражаются и заболевают, болезнь продолжается течение в течение времени $\tau$, после чего инфицированные выздоравливают, приобретая иммунитет, или гибнут (Recovered/Removed), согласно схеме:
Таким образом, кинетическая схема закона закона действующих масс имеет вид $$S+I\stackrel{k}{\rightarrow} 2I,$$ $$I\stackrel{\tau^{-1}}{\rightarrow}R.$$
Стандартный метод перевода кинетической схемы в систему обыкновенных дифференциальных уравнений приводит к
\begin{eqnarray*}\frac{dS}{dt}&=&-kSI,\\
\frac{dI}{dt}&=&kSI-\tau^{-1}I,\\
\frac{dR}{dt}&=&\tau^{-1}I.
\label{SIR}
\end{eqnarray*}
Достаточно подробное и популярное описание этой системы, а также онлайн-симуляторы ее решения могут быть найдены по ссылке.
Важно, что уравнения SIR-модели могут быть проанализированы аналитически, сводясь к одному нелинейному дифференциальному уравнению, которое, в приближении слабой эпидемии приближенно приводится к уже известному и изученному уравнению Ферхюльста для популяционного роста.
Такая редукция основана на том, что SIR-система имеет два закона сохранения (что означает алгебраическую связь между переменными, т.е. не все три из них независимы).
Первый закон сохранения - сохранение полной численности популяции (здесь $S$, $I$, $R$ - доли от полного числа, принятого за единицу, т.е. эпидемический процесс просто переводит индивидов между группам): $S+I+R=1$. Справедливость этого равенства проверятся просто - достаточно сложить три уравнения SIR модели, получится $\frac{d}{dt}(S+I+R)=0$, производная константы равно нулю, т.е. $S+I+R=\mathrm{const}$, и в качестве константы берется единица, если выражать отдельные переменные как доли полной популяции.
Второй закон сохранения получатся, если исключить квадратичный нелинейный член, "разделив" уравнения для $S$ на уравнение для $R$ ($I$ при этом сократится):
$$\frac{dS}{dR}=-k\tau S.$$
Далее, методом разделения переменных
$$\frac{dS}{S}=-k\tau dR,$$
что после интегрирования дает $\log(S)+k\tau R=\log(S_0)$, где $S_0$ - начальная доля восприимчивых особей (начальное $R_0=0$ так как в самый первый момент времени еще никто не успел переболеть).
Выразив $S$ через $R$ и подставив в $I=1-S-R$, используемое как правая часть для исходного уравнения для $dR/dt$, в результате получаем одно нелинейное уравнение
$$\frac{dR}{dt}=\tau^{-1}\left[1-R-S_0e^{-k\tau R}\right].
$$
При условии слабой эпидемии (то есть, когда переболевает далеко не вся популяция $k\tau \mathrm{max}(R)<<1$, изначально состоящая из полностью восприимчивых индивидов $S_0\approx 1$), экспоненциальный член можно разложить в ряд Тейлора до квадратичного члена
$$
\exp(-\Re R)=1-\Re_0 R+\frac{(\Re_0 R)^2}{2}
$$
и основное уравнение упрощается до уравнения Ферхюльста
$$
\frac{dR}{dt}=\frac{\Re_0-1}{\tau}R\left[1-\frac{\Re_0^2}{2(\Re_0-1)} R\right],
$$
свойства которого хорошо известны.
В частности, из его формы видно, что коэффициент роста $r=(\Re_0-1)/\tau$ и эпидемия случается, когда базовое число воспроизведения $\Re_0>1$, а выход на стационар по числу переболевших соответствует их числу (аналог емкости среды) $R_{max}\equiv K=2(\Re_0-1)/\Re_0^2$, также зависящему от базового числа воспроизведения (т.е. интересно, что ожидаемое максимальное число переболевших можно оценить, зная коэффициент роста уже в начальной экспоненциальной фазе роста).
Имея же решение для $R$, динамику числа инфицированных можно найти простым дифференцированием $I=\tau \frac{dR}{dt}$ решения уравнения Ферхюльста (или же попросту умножив правую часть уравнения Ферхюльста на характерное время протекания болезни).
Про базовое число воспроизведения:
\exp(-\Re R)=1-\Re_0 R+\frac{(\Re_0 R)^2}{2}
$$
и основное уравнение упрощается до уравнения Ферхюльста
$$
\frac{dR}{dt}=\frac{\Re_0-1}{\tau}R\left[1-\frac{\Re_0^2}{2(\Re_0-1)} R\right],
$$
свойства которого хорошо известны.
В частности, из его формы видно, что коэффициент роста $r=(\Re_0-1)/\tau$ и эпидемия случается, когда базовое число воспроизведения $\Re_0>1$, а выход на стационар по числу переболевших соответствует их числу (аналог емкости среды) $R_{max}\equiv K=2(\Re_0-1)/\Re_0^2$, также зависящему от базового числа воспроизведения (т.е. интересно, что ожидаемое максимальное число переболевших можно оценить, зная коэффициент роста уже в начальной экспоненциальной фазе роста).
Имея же решение для $R$, динамику числа инфицированных можно найти простым дифференцированием $I=\tau \frac{dR}{dt}$ решения уравнения Ферхюльста (или же попросту умножив правую часть уравнения Ферхюльста на характерное время протекания болезни).
Рекомендуемая литература по теме занятия.
Помимо приведенной выше ссылки на онлайн-материал, SIR-модель подробно описана в книге Мюррея (см. список рекомендованных учебных материалов).
Общая схема построения SIR-модели и ее анализ:
Общая схема построения SIR-модели и ее анализ:
- Blackwood J.C., Childs L.M. An introduction to compartmental modeling for the budding infectious disease modeler. Letters in Biomathematics 5 (2018), 195-221 [ссылка на статью в свободном доступе]
- Roberts M.G. (2007). The pluses and minuses of R0. Journal of the Royal Society Interface. 4 (2007) 949-96 (ссылка на статью в свободном доступе).
- van den Driessche, P. (2017). Reproduction numbers of infectious disease models. Infectious Disease Modelling 2 (2017) 288-303 (ссылка на статью в свободном доступе).
Про метод сведения SIR модели к уравнению Ферхюльста и его практическое применение:
- Lavrova A.I., Postnikov E.B., Manicheva O.A., Vishnevsky B.I. (2017). Bi-logistic model for disease dynamics caused by Mycobacterium tuberculosis in Russia. Royal Society open science, 4 (2017) 171033 (ссылка на статью в свободном доступе).
