Пандемия COVID-19 глазами математика, или почему классическая модель SEIRD не работает

Пандемія COVID-19 очима математика, або чому класична модель SEIRD не працює

За

casbt1osint.blogspot.com

13 хв

Переглянути оригінал

Анотація, або про дозвілля молодих вчених

Останні кілька тижнів ми з колегами закінчуємо робочий день тим, що змагаємося в точності прогнозу розвитку епідемії COVID-19 в Росії, використовуючи різні методи нелінійної регресії. І якщо прогноз на завтрашній день неминуче виявляється хороший, то прогноз на термін більше одного тижня відображає реальність лише в загальних рисах. Здавалося б, все зрозуміло: є епідеміологічні моделі , є методи оптимізації , є досить докладні дані , - досить поєднати це воєдино і отримати точний прогноз на місяць, а то й півроку, вперед. У цій статті я поділюся своїми міркуваннями, що не так з класичною моделлю SEIRD і як це виправити. І, звичайно, відкрию завісу таємниці, що огортає наше з вами майбутнє.

На малюнку вище наведено загальне число підтверджених випадків COVID-19 в логарифмічному масштабі для Росії і трьох європейських країн, що входять в топ-5 за кількістю заражених. Пояснення далі в тексті.

Хвилинка турботи від НЛО

У світі офіційно оголошена пандемія COVID-19 - потенційно важкої гострої респіраторної інфекції, спричиненої коронавірусів SARS-CoV-2 (2019 nCoV). На Хабре багато інформації по цій темі - завжди пам'ятайте про те, що вона може бути як достовірною / корисної, так і навпаки.

Ми закликаємо вас критично ставитися до будь-якої інформації, що публікується інформації

Офіційні джерела

Якщо ви проживаєте не в Росії, зверніться до аналогічних сайтів вашої країни.

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

Читати публікації про: коронавірус | віддалену роботу

Модель SEIRD

Модель епідемії SEIRD відноситься до класу так званих компартментальних моделей , суть яких полягає в тому, щоб розділити популяцію на кілька груп ( англ. compartments), в нашому випадку: $ inline $ S $ inline $ ( англ. susceptible) - сприйнятливі, $ inline $ E $ inline $ ( англ. exposed) - ті, у кого хвороба знаходиться в інкубаційному періоді, $ inline $ I $ inline $ ( англ. infectious) - хворі, $ inline $ R $ inline $ ( англ. recovered) - видужали, $ inline $ D $ inline $ ( англ.dead) - померлі. Потім, чисельність кожної з груп зіставляється зі змінною в системі диференціальних рівнянь, вирішуючи яку, можна спрогнозувати динаміку розвитку епідемії. Модифікацій моделі SEIRD досить багато, наприклад, SEIR - спрощена модель, яка не враховує окремо випадки одужання і смерті. Для ознайомлення з іншими моделями можу порекомендувати непогану статтю на цю тему.

трохи теорії

Вперше модель епідемії у вигляді системи з трьох диференціальних рівнянь для змінних $ inline $ S, I, R $ inline $ з'явилася в роботі У. Кермак і А. Мак-Кендрика 1927 року.

Ці диференціальні рівняння мають вигляд:

$$ відображення $$ \ початок {вирівнювання} \ frac {dS} {dt} & = - \ beta \ frac {SI} {N}, \\ \ frac {dI} {dt} & = \ beta \ frac {SI } {N} - \ gamma I, \\ \ frac {dR} {dt} & = \ gamma I, \ end {align} $$ display $$

де, крім знайомих нам змінних фігурують такі константи: $ inline $ N = S + I + R $ inline $ - загальний розмір популяції, $ inline $ \ beta $ inline $ - швидкість передачі інфекції, $ inline $ \ gamma $ inline $ - швидкість одужання.

Сенс рівняння Кермак і Мак-Кендрика наступний: число сприйнятливих убуває пропорційно їх числа, помноженому на середню частку інфікованих в популяції $ inline $ I / N $ inline $, число інфікованих приростає тими ж темпами з поправкою на те, що деяке їх число $ inline $ \ gamma I $ inline $ одужує, і число видужали приростає за рахунок зменшення числа інфікованих. Варто відзначити, що модель SIR містить нелінійність $ inline $ SI $ inline $, через що аналітичне рішення системи рівнянь стає в загальному випадку неможливим, але, благо, методи чисельного диференціювання легко справляються з цим завданням.

Додавши сюди ще одну змінну $ inline $ E $ inline $ (число людей з хворобою в інкубаційному періоді), отримаємо SEIR-модель:

$$ відображення $$ \ початок {вирівняти} \ frac {dS} {dt} & = - \ beta \ frac {SI} {N}, \\ \ frac {dE} {dt} & = \ beta \ frac {SI } {N} - \ kappa E, \\ \ frac {dI} {dt} & = \ kappa E- \ gamma I, \\ \ frac {dR} {dt} & = \ gamma I, \ end {align} $$ показ $$

де з'являється ще одна константа $ inline $ \ kappa $ inline $ - швидкість переходу хвороби з інкубаційної стадії в відкриту. Малюнок з узятий з статті .

Відразу видно, що SEIR-модель не дуже годиться для опису COVID-19 хоча б тому, що в цій моделі приховані носії інфекції $ inline $ E $ inline $ незаразних. Виправити цей недолік можна, запровадивши, слідом за Pengpeng et al. , Додатковий параметр $ inline $ \ theta $ inline $, що характеризує ступінь заразність латентних носіїв інфекції в порівнянні з хворими. Модифікована модель SEIR, яку ми спробуємо застосувати до поточної епідемії, матиме вигляд:

$$ display $$ \ begin {align} \ frac {dS} {dt} & = - \ beta \ frac {S (I + \ theta E)} {N}, \\ \ frac {dE} {dt} & = \ beta \ frac {S (I + \ theta E)} {N} - \ kappa E, \\ \ frac {dI} {dt} & = \ kappa E- \ gamma I, \\ \ frac {dR} {dt} & = \ gamma I. \ end {align} $$ display $$

На перший погляд, отримана модель обіцяє бути цілком правдоподібною.

Чисельний експеримент з моделлю SEIR

Для моделювання спробуємо взяти такі параметри, орієнтуючись на відкриті дані . Припускаючи, що хвороба в середньому триває 14 днів (по крайней мере, скільки триває легка форма, на яку припадає до 80% випадків), знайдемо значення $ inline $ \ gamma = 1/14 = 0,0714 $ inline $. Приймемо $ inline $ \ beta = 3/14 = 0,2143 $ inline $. Величину $ inline $ \ theta = 0,6 $ inline $ запозичуємо з роботи Pengpeng et al . C урахуванням середньої тривалості інкубаційного періоду в 3 дня, візьмемо $ inline $ \ kappa = 1/3 = 0,33 $ inline $. Населення Росії приймемо рівним $ inline $ N = 144,5 \ cdot10 ^ 6 $ inline $ людина.

В якості вихідних умов використовуємо дані по Росії на 2 квітня, коли введені в кінці березня заходи щодо обмеження поширення інфекції повинні були здобути свою дію, а саме:

$$ відображення $$ \ begin {align} & S_0 = 3548, \\ & I_0 = 3283, \\ & E_0 = 0,5 I_0. \ end {align} $$ показ $$

Оцінку $ inline $ E_0 $ inline $ ми взяли щодо довільно, та це й неважливо, оскільки, як ви розумієте, що щось піде не так.

В результаті моделювання методом Ейлера з кроком в 1 день з 2-го по 24 квітня включно, одержимо графіки, наведені нижче: зліва в лінійному масштабі, праворуч в логарифмічному.

Круглими маркерами відзначені реальні дані за загальною кількістю випадків в Росії, квадратними - по числу хворих. На перший погляд, результати виглядають непогано, крім одного: з параметрами моделі ми явно не вгадали. І тут нам на допомогу приходять методи оптимізації.

Оптимізуй це

Методи оптимізації, якщо читач з ними не знайомий, - це алгоритми, що дозволяють відшукати мінімум деякої цільової функції. У нашому випадку перед нами - завдання нелінійної регресії: як підібрати вектор параметрів диференціального рівняння $ inline $ \ mathbf {x} = (\ beta, \ gamma, \ kappa, E_0) ^ \ top $ inline $ так, щоб набір точок рішення диференціального рівняння $ inline $ F $ inline $ був максимально близький набору точок спостереження $ inline $ \ mathcal {F} $ inline $.

Скористаємося середньоквадратичним відхиленням як мірою похибки моделі. Цільова функція набуде вигляду

$$ відображення $$ f (\ mathbf {x}) = \ frac {1} {M} \ sqrt {\ sum_ {i = 1} ^ {M} (F_ {i} - \ mathcal {F} _ {i }) ^ 2 + \ sum_ {i = 1} ^ {M} (G_ {i} - \ mathcal {G} _ {i}) ^ 2}, $$ відображення $$

де $ inline $ M $ inline $ - число точок, $ inline $ F $ inline $ - загальне число випадків зараження, яке дає модель, $ inline $ \ mathcal {F} $ inline $ - реальне загальне число випадків, $ inline $ G $ inline $ - число хворих в поточний момент, який дає модель, $ inline $ \ mathcal {G} $ inline $ - реальне загальне поточних активних випадків.

Скориставшись Optimization Toolbox в MATLAB, підгонимо параметри моделі під дані спостереження. В результаті отримаємо рішення, наведене на малюнку нижче.

На перший погляд, все відмінно. Невязка вийшла рівною $ inline $ f (\ mathbf {x}) = 131,98 $ inline $, та й «на опуклий морське око» підгонка рішення вийшла на славу. Подивимося на отримані параметри:

$$ display $$ \ begin {align} & \ beta = 0,374, \\ & \ gamma = 0,0117, \\ & E_0 = 7,84 \ cdot 10 ^ 6, \\ & \ kappa = 4,81 \ cdot 10 ^ {- 5}. \ end {align} $$ показ $$

Величина майже 8 млн. Латентних хворих при зареєстрованих близько 60 тис. Випадків на

24 серпня - щось сумнівне. У нас також вийшло, що середній час переходу в активну фазу хвороби дорівнює $ inline $ 1 / \ kappa = 2079 $ inline $ днів.

Чому так вийшло? Все стане зрозуміло, якщо ми проаналізуємо форму кривої на тривалому масштабі часу. Для цього візьмемо нашу SEIR-модель з «правдоподібними» параметрами і Промоделюємо на тривалому проміжку часу (в цьому досвіді я прийняв нове значення $ inline $ \ beta = 0,186 $ inline $):

Крива, що відповідає загальній кількості випадків, має характерну S-образну форму в лінійному масштабі. Цю-то форму і спробувала надати кривої оптимізаційна програма. Крім того, що сам прогноз з «правдоподібними» параметрами жахає - по ньому, до вересня перехворіє майже 90% населення країни - він очевидно нереалістичний, якщо подивитися результати по іншим країнам (та ж картинка, що і на початку статті, тільки в лінійному масштабі ):

Тут я порівнюю три європейські країни, що знаходяться в топ-5 за кількістю хворих, і Росію. Видно, що за темпами розвитку епідемії ми відстаємо приблизно на місяць, і що ось уже як місяць зростання загального числа випадків у всіх трьох країнах практично лінійний (і навіть повільніше лінійного), на відміну від результатів, отриманих в SEIR моделі. Звідси виникають три питання:

1. Чому зростання епідемії уповільнюється до лінійного?
2. Як змінити класичну модель SEIR, щоб вона знову була релевантною?
3. Чому, якщо зростання епідемії лінійний, ми все одно не можемо нічого впевнено передбачити на місяць або рік вперед?

Почну з відповіді на третє питання. Коли ми щось прогнозуємо, перед нами постає завдання досить неприємна: дані, на яких ми будуємо модель, неідеальні - вони містять помилки, шум, і побудована на їх основі модель також містить деяку помилку. Коли ми продовжуємо тимчасової ряд нашими модельними точками, помилка накопичується - і досить швидко, якщо ми прогнозуємо зростаючу за часом функцію. А це має місце як раз в нашому випадку. Більш того, модель на те і модель, що відображає реальну ситуацію дуже обмежено. Раптовий розвиток епідемії в новому великому місті, застосування більш ефективного способу лікування, зміна способу збору інформації - все це може внести в реальні дані стільки помилок, що довгостроковий прогноз виявиться абсолютно далекий від реальності.

Милиці та велосипеди: модифікуємо модель SEIR

Спробуємо відповісти на питання, чому зростання епідемії уповільнюється до лінійного. При тій кількості заражених, яке ми маємо зараз, значну роль починає грати ефект масштабу, пов'язаний з обмеженою швидкістю комунікацій між людьми.

Якщо говорити більш точно, то згадаємо: кількість хворих в моделі SEIR прямо пропорційно середньому числу хворих в популяції $ inline $ I / N $ inline $. Це правило добре працює в невеликих популяціях, де кожен може спілкуватися з кожним, а хворі розподілені рівномірно. У реальності, особливо в масштабах десятків і сотень тисяч людей, якщо взяти двох випадкових хворих людей, то виявиться, що вони не тільки ніколи один з одним не спілкувалися і не бачили один одного, вони навіть не їздили в одному і тому ж вагоні метро. Та й взагалі живуть в різних містах. Все, що їх об'єднує - ланцюжок соціальних зв'язків, яка призвела до того, що їм передався вірус.

Як приклад я побудував модель епідемії у вигляді клітинного автомата, де кожна клітина взаємодіє тільки з 4 сусідніми. Це еквівалентно тому, що у кожного індивіда популяції 4 соціальних контакту - це дуже маленьке число для людської популяції, але тим швидше проявиться ефект обмеження соціальних зв'язків. На кожній ітерації з ймовірністю 0,1 кожен з 4-х сусідів зараженої клітини може бути заражений. Хвороба триває в середньому 14 днів. Результати моделювання для пулу з 200x200 клітин представлені на малюнку нижче, де $ inline $ k $ inline $ - номер ітерації.

Синій колір позначає сприйнятливих, жовтий - хворих, зелений - видужали. Найцікавіше - як виглядають графіки чисельності хворих. А виглядають вони приблизно так, як і планувалося: після недовгої фази субекспоненціальное зростання, прямо як в моделі SEIR, спостерігається затяжна фаза лінійного зростання - прямо як в реальності.

У мене не було мети отримати картинку, схожу на реальність кількісно. Якщо хочеться більшої правдоподібності - можу порекомендувати проект Сергія Потєхіна, про який була недавно публікація на Хабре . Для допитливих читачів нижче наведено більш суворе доказ лінійності зростання.

Доказ теореми про лінійному зростанні епідемії в великій популяції

Візьмемо геометричну інтерпретацію: нехай граф соціальних зв'язків представлений у вигляді $ inline $ d $ inline $ -мірною решітки. У моделі у вигляді клітинного автомата решітка двовимірна. У реальності, при середньому числі щоденних соціальних контактів в 20 (оцінка з вищезгаданої публікації ) розмірність може бути грубо оцінена як $ inline $ d \ approx 4 $ inline $. Кожен носій інфекції породжує зростаючий навколо нього $ inline $ d $ inline $ -мірною гиперкуб з вдруге заражених. Ребро куба має довжину $ inline $ n ^ {\ frac {1} {d}} $ inline $, і, якщо контакт із зараженою призводить до захворювання з імовірністю $ inline $ P $ inline $, то кожен день кожне ребро в середньому подовжується на $ inline $ 2P $ inline $. Таким чином, отримуємо модель зростання, що виражається рекурентним співвідношенням

$$ відображення $$ n_ {k + 1} = \ ліворуч (n_k ^ {\ frac {1} {d}} + 2P \ праворуч) ^ d, $$ відображення $$

Уявімо $ inline $ n $ inline $ як функцію часу: $ inline $ n_k = n (t_k); n_ {k + 1} = n (t_k + 1) $ inline $. Продифференцировав попереднє співвідношення, отримаємо:

$$ відображення $$ n '_ {k + 1} = n' _ {k} \ зліва (1 + \ frac {2P} {n_k ^ {\ frac {1} {d}}} \ праворуч) {{d -1}. $$ відображення $$

Звідси видно, що при великих значеннях $ inline $ n_k $ inline $ похідні дорівнюють $ inline $ n '_ {k + 1} = n' _ {k} $ inline $, отже, зростання лінійний.

Нижче представлено розвиток сценарію при ймовірності заразитися 4% і 16-ти соціальних зв'язках одного індивіда.

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

Завдання для допитливих

Тепер про модифікацію моделі SEIR. Найпростіше, що ми можемо зробити - помножити нелінійний компонент на деяку функцію, залежну від кількості хворих. При малому числі хворих ця функція повинна бути близька до 1, при великому - повинна асимптотично прагне до нуля. Найпростішим підходящим кандидатом є

$$ відображення $$ \ mathcal {\ varphi} (I, E) = e ^ {- \ alpha (I + \ theta E) ^ {K_0}}. $$ display $$

Підбором параметрів $ inline $ \ alpha $ inline $ і $ inline $ K_0 $ inline $ можна компенсувати експоненціальне зростання в оригінальної моделі.

Додамо в модель, для більшої інформативності, і компоненту $ inline $ D $ inline $ - число смертей. Отримаємо модифікований варіант SEIRD-моделі:

$$ display $$ \ begin {align} \ frac {dS} {dt} & = - \ beta \ frac {S (I + \ theta E) \ mathcal {\ varphi} (I, E)} {N}, \\ \ frac {dE} {dt} & = \ beta \ frac {S (I + \ theta E) \ mathcal {\ varphi} (I, E)} {N} - \ kappa E, \\ \ frac { dI} {dt} & = \ kappa E- \ gamma I- \ mu D, \\ \ frac {dR} {dt} & = \ gamma I, \\ \ frac {dD} {dt} & = \ mu D . \ end {align} $$ показ $$

Результати моделювання показані на малюнку нижче.

Ефективне значення похибка в порівнянні з оригінальною моделлю майже не змінилася. Величини параметрів вже реалістичні. Для зручності я позначив початкове число заражених в активній фазі як $ inline $ I_0 $ inline $.

$$ display $$ \ begin {align} & \ beta = 0,219, \\ & \ gamma = 0,0102, \\ & E_0 = 0,13 \ cdot I_0, \\ & \ kappa = 1/3, \\ & \ mu = 1,13 \ cdot 10 ^ {- 3}. \ end {align} $$ показ $$

Модель дуже добре интерполирует і похідні - величини щоденного приросту числа хворих і кількості смертей.

Спробуємо зробити прогноз. Візьмемо горизонт прогнозування 2 місяці і продовжимо рішення моделі зі знайденими оптимизационной програмою параметрами.

На перший погляд, непогано, але такий прогноз улюбленій батьківщині не побажаєш: число нових випадків буде продовжувати знижуватися, але загальне число продовжить рости. Зупинити епідемію в цьому випадку можна тільки за допомогою вакцини, або почекавши, поки перехворіє майже все населення. Число нових смертей встановлюється рівним приблизно 200 в день. Це наочна ілюстрація того, що буде, якщо не посилювати заходи по боротьбі з епідемією. Невже нас чекає саме це? І заради цього не дуже світлого майбутнього багато з нас старанно сидять вдома, закупившись гречкою і туалетним папером?

Нижче я розгляну два сценарії, і дивлячись на туманну далечінь прийдешніх місяців із 28 квітня 2020 року, не можу точно сказати, яким із них будуть розвиватися події далі. Зараз, в момент перелому кривої нових випадків, ми знаходимося в точці, звідки що-небудь передбачити подвійно проблематично.

Сценарий США

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

Модифікована модель SEIRD, навчена на перших 33 точках, починаючи з 2 березня, плюс-мінус реалістично пророкує протягом епідемії в квітні.

Як бачимо, зростання в квітні практично ідеально лінійний. Модель трохи завищує смертність за квітень, але загальна картина виявляється вірною.

На цій картинці показані дані по щоденному приросту нових випадків і смертей в США. Дуже схоже на те, що передбачала модель для Росії.

сценарій Німеччини

Дисципліновані німці зуміли переломити хід кривої в свою користь, і її зростання відбувається повільніше лінійного. Більш того, щоб зробити модель релевантної, мені довелося вручну додати 6 квітня збільшення коефіцієнта одужання $ inline $ \ gamma $ inline $ в 1,7 раз, інакше таке різке падіння числа випадків в термінах моделі SEIRD не пояснити.

Модель навчалася на перших 27 точках, починаючи з 10 березня. Також я змінив і нелінійну функцію. Для Німеччини краще підійшла експонента, що залежить від часу:

$$ display $$ \ mathcal {\ varphi} (t) = K_0 e ^ {- \ alpha t}. $$ display $$

Такий вид функції свідчить про кумулятивний наростанні числа перерваних соціальних зв'язків і, відповідно, шляхів поширення інфекції. Ось вам і наочна ілюстрація користі самоізоляції.

Вище показані величини щоденного приросту нових випадків і смертей. Як і в випадку США, реальні дані містять ясно виражені коливання з періодом в 7 днів. Це означає, що у вихідні дні число контактів збільшується, а отже, зростає і число заражених.

висновок

Робити прогнози - як короткострокові, так і довгострокові - не просто данина цікавості. У разі епідемії потрібно знати, скільки ліжкомісць слід підготувати, скільки апаратів ШВЛ зробити, на скільки місяців зробити запаси засобів індивідуального захисту для медиків. Посадові особи повинні розуміти, чи достатньо вжитих заходів, або слід ввести нові заборони і обмеження. В ідеалі, модель повинна настільки добре відображати реальність, щоб по ній можна було побачити силу дії кожної знову прийнятої заходи, і тоді можна було б посилювати корисні заходи і скасовувати рішення, які опинилися марними.

Хоч і з деякими застереженнями, але всіх, хто поки ще не набув імунітет від COVID-19, дійсно можна описати буквою $ inline $ S $ inline $, все різноманіття хворих людей - буквою $ inline $ I $ inline $, і так далі. Більш того, за допомогою моделі SEIRD можна навіть дещо пояснити. Але передбачити щось у віддаленому майбутньому вона може вкрай приблизно.
Я навмисне навів у статті тільки негативний сценарій, коли ми повторюємо долю США в плані динаміки епідемії. Якщо цей сценарій виявиться вірним, до кінця червня у нас буде більше 300 тис. Зареєстрованих випадків захворювання і більше 10 тисяч смертей. Хоча є передумови до того, що цей сценарій не втілиться в реальність, я б порадив поставитися до нього за принципом: «сподівайся на краще, готуйся до гіршого». Як то кажуть, якщо вже за боротьбу з епідемією в США взялися кращі уми НАСА , значить, справа і справді погане.
Поки все, що нам залишається - менше відвідувати громадські місця, користуватися правильними респіраторами , мити руки, не забувати протирати спиртом смартфон після того, як дістали його на вулиці, і дотримуватися інші прості рекомендації.

Але все ж, чи є можливість робити більш точні прогнози? Так звичайно. Але про це якось іншим разом.

Якщо у вас є бажання пограти з вихідним кодом і самим запропонувати варіант розвитку подій, то ось посилання на гітхаб .

Просмотры: