Моделирование эпидемии: модель SIR
Как тремя уравнениями описать распространение болезни, предсказать пик заражений и понять, почему число R0 решает всё.
Модель SIR — компартментная модель эпидемии, делящая население на восприимчивых (S), инфицированных (I) и выздоровевших/удалённых (R), и описывающая переходы между группами дифференциальными уравнениями.
Во время любой эпидемии главный вопрос — насколько быстро болезнь распространится, когда наступит пик и сколько людей переболеет. Удивительно, но грубый, но очень информативный ответ даёт модель из трёх уравнений, предложенная Кермаком и Маккендриком ещё в 1927 году. Она лежит в основе того, что эпидемиологи рассказывали по телевизору во время пандемий.
Уравнения модели SIR
Население N делится на три группы. S — susceptible, восприимчивые (могут заразиться). I — infected, заражённые (распространяют болезнь). R — recovered, выздоровевшие или умершие (больше не участвуют). Уравнения:
S' = -beta * S * I / N I' = beta * S * I / N - gamma * I R' = gamma * I
Параметр beta — скорость передачи (сколько контактов с заражением в единицу времени), gamma — скорость выздоровления (доля инфицированных, выздоравливающих за единицу времени). Член beta*S*I/N — это число новых заражений: оно пропорционально и числу восприимчивых, и числу инфицированных. Член gamma*I — число выздоровлений. Заметьте: S+I+R всегда равно N, потому что S'+I'+R'=0. Никто не теряется и не появляется.
Базовое репродуктивное число R0
Ключевая величина модели — R0 = beta / gamma. Это среднее число людей, которых заразит один больной в полностью восприимчивом окружении. Смысл порога прост: в начале эпидемии S примерно равно N, поэтому I' приблизительно равно (beta - gamma)*I. Если beta > gamma, то есть R0 > 1, число инфицированных растёт — начинается эпидемия. Если R0 < 1, болезнь затухает сразу. Именно поэтому вакцинация и карантин нацелены на снижение R0 ниже единицы: тогда вспышка гаснет сама.
Из R0 вытекает понятие порога коллективного иммунитета — той доли населения, которую достаточно сделать невосприимчивой (переболевшей или привитой), чтобы эпидемия не могла разрастаться. Этот порог равен 1 − 1/R0. Для нашего R0 = 3 он составляет 1 − 1/3, то есть примерно 0.67 — около 67% популяции. Логика прямая: как только две трети окружения больного уже невосприимчивы, в среднем он передаёт болезнь меньше чем одному новому человеку, и цепочка заражений рвётся. Чем заразнее болезнь и чем больше R0, тем выше порог и тем больше людей должно получить иммунитет, чтобы защитить и тех, кто привиться не может.
Важный и не сразу очевидный момент: по ходу эпидемии «работающее» репродуктивное число — его называют эффективным, Rt — не остаётся постоянным, а падает. Каждый переболевший уходит из группы S, и доля восприимчивых вокруг любого больного сокращается. Фактически Rt = R0 · S/N: пока восприимчивых много, болезнь летит, но по мере «выгорания» группы S множитель S/N снижается, и в какой-то момент Rt опускается до единицы. Именно тогда — на пике — рост сменяется спадом. Поэтому эпидемия останавливается до того, как переболеют все: ей просто перестаёт хватать восприимчивых, чтобы поддерживать саму себя.
Как работает под капотом
Численно мы решаем систему так же, как Лоренца: состояние [S, I, R], правая часть возвращает три производные, RK4 продвигает на шаг h. Возьмём население N=1000, начнём с одного больного (I0=1, S0=999, R0_start=0), beta=0.3, gamma=0.1, то есть R0 = 0.3/0.1 = 3. Это значит, что без мер каждый больной заражает в среднем троих — эпидемия неизбежна. Прогоним 80 дней с шагом h=0.5 и распечатаем состояние каждые 10 дней.
def sir(state, beta=0.3, gamma=0.1, N=1000.0):
S, I, R = state
new_inf = beta * S * I / N
rec = gamma * I
return [-new_inf, new_inf - rec, rec]
def rk4_step(f, state, h):
k1 = f(state)
s2 = [state[i] + 0.5 * h * k1[i] for i in range(3)]
k2 = f(s2)
s3 = [state[i] + 0.5 * h * k2[i] for i in range(3)]
k3 = f(s3)
s4 = [state[i] + h * k3[i] for i in range(3)]
k4 = f(s4)
return [state[i] + (h/6.0)*(k1[i] + 2*k2[i] + 2*k3[i] + k4[i]) for i in range(3)]
state = [999.0, 1.0, 0.0]
h = 0.5
print("день S I R")
print(str(0).rjust(3), " ", round(state[0],1), " ", round(state[1],1), " ", round(state[2],1))
day = 0.0
for step in range(1, 161):
state = rk4_step(sir, state, h)
day += h
if abs(day - round(day)) < 1e-9 and int(round(day)) % 10 == 0:
d = int(round(day))
print(str(d).rjust(3), " ", round(state[0],1), " ", round(state[1],1), " ", round(state[2],1))Вывод:
день S I R 0 999.0 1.0 0.0 10 963.5 27.3 9.2 20 642.4 222.8 134.8 30 201.1 353.6 445.3 40 91.3 180.1 728.6 50 67.5 74.0 858.5 60 59.5 28.6 911.9 70 56.2 10.8 933.0 80 54.9 4.0 941.1
Видна классическая картина: число инфицированных I сначала растёт, достигает пика где-то около 30-го дня (порядка 350 человек одновременно болеют), затем спадает. К концу почти все оказываются в группе R (переболели), и остаётся около 55 человек, которые так и не заразились. Это и есть «коллективный иммунитет»: эпидемия гаснет не потому, что закончились больные, а потому, что закончились восприимчивые.
ASCII-кривая I(t)
Нарисуем грубую кривую числа инфицированных по дням, чтобы увидеть пик глазами. Каждая звёздочка — примерно 8 человек:
день | I(t) 0 | 10 | *** 20 | **************************** 30 | ******************************************** 40 | ********************** 50 | ********* 60 | *** 70 | * 80 |
Форма «горба» — фирменный знак эпидемии. Цель всех мер сдерживания (масок, карантина, дистанции) — снизить beta, тем самым «сгладить кривую»: сделать пик ниже и шире, чтобы больницы справились с нагрузкой.
Частые ошибки
- Забытое деление на N. Член заражения — это beta*S*I/N, а не beta*S*I. Без деления масштаб скорости получится в N раз больше и эпидемия «вспыхнет» мгновенно.
- I0 = 0. Если стартовать без единого больного, I' = 0 навсегда — эпидемии не будет. Нужен хотя бы один заражённый.
- Путаница R (выздоровевшие) и R0 (репродуктивное число). Это совершенно разные вещи с похожим обозначением.
- Ожидание, что переболеют все. Даже при R0=3 часть населения (здесь ~55 человек) никогда не заражается — это нормальное свойство модели.
Зачем это нужно
SIR — фундамент математической эпидемиологии. Её расширения (SEIR с инкубационным периодом, модели с возрастными группами и вакцинацией) использовались для планирования мер во время реальных пандемий. Понимание R0 и порога эпидемии объясняет, почему важен охват вакцинацией и как работает коллективный иммунитет. А с точки зрения численных методов это идеальный пример системы ОДУ, где маленькие уравнения дают богатую и практически важную динамику.
Полезно ясно понимать и границы применимости базовой модели — именно из её упрощений и рождаются все расширения. SIR предполагает однородное перемешивание: считается, что любой больной с равной вероятностью контактирует с любым человеком в популяции, как молекулы в хорошо размешанном растворе. В реальности контакты структурированы — семья, школа, работа, город, — и болезнь распространяется по сети связей неравномерно. Модель также не различает возрастные группы, хотя для многих инфекций риск и заразность сильно зависят от возраста. Наконец, в SIR человек заразен сразу после контакта, тогда как у реальных болезней есть латентный (инкубационный) период. Чтобы учесть его, добавляют четвёртую группу E (exposed, заражён, но ещё не заразен) — так получается модель SEIR; дальше идут модели с возрастными контактными матрицами, с пространственной структурой и с повторным заражением.
Через ту же модель удобно рассуждать о мерах. Любое вмешательство в конечном счёте меняет два числа — beta или gamma. Карантин, маски, дистанция и закрытие массовых мероприятий снижают beta, потому что урезают число «контактов с заражением» в единицу времени; именно так пик «сглаживается». Быстрое тестирование и изоляция заболевших, по сути, увеличивают gamma — больной раньше выводится из игры и меньше успевает заразить. Вакцинация действует иначе: она не трогает beta и gamma напрямую, а сразу переводит часть людей из S в R, ещё до болезни, уменьшая долю восприимчивых и тем самым снижая эффективное Rt. Все эти рычаги бьют в одну цель — увести систему за порог, где R0 (или хотя бы Rt) опускается ниже единицы.
- SIR делит население на восприимчивых, инфицированных и выздоровевших; S+I+R = N сохраняется.
- R0 = beta/gamma — порог: при R0 > 1 эпидемия растёт, при R0 < 1 затухает.
- Численно решается RK4 для системы [S, I, R]; кривая I(t) имеет характерный пик.
- Меры сдерживания снижают beta и 'сглаживают кривую', уменьшая пиковую нагрузку.