Pomiar i prognoza pandemii COVID-19 w Polsce w czasie rzeczywistym
Rafał J. Mostowy
Małopolskie Centrum Biotechnologii Uniwersytetu Jagiellońskiego, Kraków, Polska
Data publikacji: 30/03/2020
Data ostatniej aktualizacji: 27/03/2020
Adres korespondencyjny: rafal.mostowy@uj.edu.pl
Podsumowanie
Polska, podobnie jak inne kraje, przechodzi obecnie pandemię COVID-19 spowodowaną wirusem SARS-CoV-2, jednak do tej pory przeprowadzono bardzo niewiele prognoz dotyczących przebiegu epidemii w Polsce. W celu pomiaru i prognozy dynamiki epidemiologicznej COVID-19 w Polsce postanowiłem wykorzystać podejście oparte na modelu matematycznym, który symuluje epidemię typu SARS i dopasowuję go danych opartych na raportowanej liczbie zgonów pacjentów COVID-19 przedstawionej przez Ministerstwo Zdrowia. Dzięki temu podejściu wyliczam prawdopodobny rzeczywisty początek epidemii w Polsce, mierzę skalę niedoszacowania liczby zakażonych i stawiam prognozę epidemiologiczną na najbliższe miesiące w zależności od skuteczności rządowych działań na rzecz dystansowania społecznego. Analizy te prowadzą do następujących ustaleń i wniosków.
Szacuję, że epidemia w Polsce rozpoczęła się mniej więcej w drugiej połowie stycznia. To sugeruje, że jest mało prawdopodobne, iż pierwszy polski odnotowany przypadek zakażenia COVID-19 w dniu 04/03/2020 jest pierwszym zakażeniem COVID-19 w Polsce, i że w tym dniu takich zakażeń było dużo więcej.
Obecnie mamy duży stopień niedoszacowania przypadków. Szacuję, że - podobnie jak w wielu innych krajach - w oficjalnych raportach brakuje nam od 50% do 75% rzeczywistych przypadków. Liczba ta pokrywa się z wcześniejszymi szacunkami dotyczącymi odsetka populacji, dla której zakażenie SARS-CoV-2 jest bezobjawowe lub łagodne.
Stopień niedoszacowania liczby zakażonych w ostatnich dniach się zmniejsza, prawdopodobnie dzięki rosnącej liczbie testów RT-PCR wykonywanych codziennie, ale tendencja w rosnącej proporcji wykrytych zakażeń będzie się utrzymywać tylko wtedy, gdy działania na rzecz dystansowania społecznego znacznie ograniczą transmisje wirusa w populacji.
Skuteczność rządowych środków dystansowania społecznego będzie miała kluczowe znaczenie dla powstrzymania epidemii w nadchodzących miesiącach, pomagając zapobiec dziesiątkom tysięcy niepotrzebnych zgonów.
Nawet jeżeli dzięki drastycznemu ograniczeniu rozprzestrzeniania się wirusa za pomocą dystansowania społecznego fala epidemiologiczna w najbliższych tygodniach zniknie, większość populacji pozostanie podatna na infekcje. Stworzy to ryzyko wybuchu kolejnej epidemii w momencie, w którym obostrzenia rządowe ograniczające kontak między ludźmi zostaną zniesione. Biorąc pod uwagę jak mało obecnie wiemy o wpływie zmian pogodowych na transmisję wirusa SARS-CoV-2, Polska potrzebuje długoterminowej strategii kontroli epidemii przez zapobieganie potencjalnym wybuchom w przyszłości do czasu wejścia na rynek ogólno-dostępnej szczepionki.
Wprowadzenie
Na dzień dzisiejszy widziałem niewiele analiz zorientowanych na predykcję dynamiki epidemiologicznej pandemii COVID-19 w Polsce. Te, które widziałem, opierają się na dopasowaniu funkcji wykładniczych do liczby zgłaszanych przypadków. Takie podejście ma co najmniej dwie poważne wady. Pierwsza wada polega na tym, że faza wykładnicza epidemii utrzymuje się tylko na początku, gdy większość populacji jest nadal podatna na zakażenie. Przez to analiz takich nie można wykorzystać do predykcji tego, kiedy krzywa epidemiologiczna zwolni i osiągnie maksimum liczby zachorowań. Drugą wadą jest to, że ponieważ przebieg choroby COVID-19 jest najprawdopodobniej za bezobjawowy w dużej części populacji, liczba zgłaszanych przypadków będzie zależała nie tylko od dynamiki epidemiologicznej, ale także od proporcji populacji, która jest testowana na obecność wirusa. Dlatego możliwość przewidywania skutków epidemii w Polsce, jak w każdym innym kraju, wymaga podejścia statystycznego, które przy absolutnym minimum jest w stanie uporać się z tymi dwoma problemami.
W tym raporcie postanowiłem posłużyć się podejściem podobnym do wypracowanego przez moich kolegów z Uniwersytetu w Bernie, które opiera się na idei dopasowania modelu epidemiologicznego typu SEIR do liczby zgłaszanych zgonów u pacjentów COVID-19, argumentując, że takie przypadki mają małe prawdopodobieństwo pominięcia. Dopasowuję ten model do oficjalnych danych Ministerstwa Zdrowia, zebranych i upublicznionych przez pana Michała Rogalskiego pod tym linkiem. Dane te wskazują, że pierwszy zgłoszony i potwierdzony przypadek zakażenia SARS-CoV-2 w Polsce miał miejsce 4 marca. Polski rząd wprowadził pierwsze ograniczenia (tzw. stan zagrożenia epidemicznego) 10 dni później, 14 marca. Korzystając z tych narzędzi, starałem się wykorzystać podejście oparte na modelowaniu matematycznym, aby uzyskać wgląd w czas początku polskiej epidemii, oszacować odsetek niewykrytych przypadków zakażenia na 26/03 oraz przewidzieć różne scenariusze epidemii w zależności od skuteczności działań ograniczeń rządowych. Na koniec otwarcie dyskutuję również ograniczenia mojego podejścia i komentuję ograniczenia modelowania matematycznego w epidemiologii. W załączeniu do raportu przedstawiam kod napisany za pomocą R MarkdownMarkdown, który każdy może pobrać i wykorzystać do dalszych analiz.
Metodologia
Analiza opiera się na matematycznym modelu epidemiologicznym typu SEIR, z istotną różnicą, że model dodatkowo uwzględnia osoby, które są hospitalizowane, przyjmowane na oddział intensywnej terapii, oraz które umierają od infekcji. Model można zwięźle podsumować następującym diagramem.
Rycina 1. Podsumowanie modelu epidemiologicznego użytego w analizje Podatni na infekcję (S) są narażeni (E) kiedy wchodzą w kontakt z osobą zakaźną (I) ze współczynnikiem \(\beta I\), a narażeni stają się zakaźni ze współczynnikiem \(\sigma\). Zakażeni zdrowieją (R) ze współczynnikiem \((1-\epsilon_1)\gamma\) i są poddawani hospitalizacji (H) ze współczynnikiem \(\epsilon_1\gamma\). Hospitalizowani zdrowieją ze współczynnikiem \((1-\epsilon_2)\omega_1\) i są przyjmowani na oddział intensywnej terapii (V) ze współczynnikiem \(\epsilon_2\omega_1\). Pacjenci oddziału intensywnej terapii zdrowieją ze współczynnikiem \((1-\epsilon_3)\omega_2\) i umierają (D) ze współczynnikiem \(\epsilon_3\omega_2\). W tym modelu \(\beta=R_0\gamma/N\), gdzie \(R_0\) to bazowy współczynnik reprodukcji a \(N\) to wielkość populacji. Epidemia zaczyna się w czasie \(t_{case}-T_{seed}\), gdzie \(t_{case}\) to data pierwszego potwierdzonego zakażenia (04/03) a \(T_\text{seed}\) to odstęp czasu między pierwszym potwierdzonym zakażeniem a pierwszym rzeczywistym zakażeniem. W dniu \(t_{case}+10\), czyli 14/03, wprowadzone są rządowe interwencje i bazowy współczynnik reprodukcji zostaje zredukowany do \(\kappa R_0\), gdzie \(\kappa\in[0,1]\).
Współczynniki w powyższym modelu zostały obliczone na podstawie założonych parametrów modelu, których wartości przyjąłem na podstawie literatury.
| Parametr | Wartość | Źródło |
|---|---|---|
| Wielkość populacji RP | 38 386 000 | Główny Urząd Statystyczny (GUS) Office |
| Interwał seryjny (okres przejściowy+zakaźny) | 7.5 dnia | Li et al. |
| Czas trwania pobytu w szpitalu w przypadku lekkich i ciężkich przypadków | 8 dni | Przyjęte na podstawie Imperial College COVID-19 Response Team: Report 9 |
| Dodatkowy czas pobytu w szpitalu w przypadkach krytycznych | 8 dni | Przyjęte na podstawie Imperial College COVID-19 Response Team: Report 9 |
| Proporcja przypadków hospitalizowanych | 5% | Przyjęte na podstawie Imperial College COVID-19 Response Team: Report 9 |
| Proporcja przypadków krytycznych | 2.5% | Przyjęte na podstawie Imperial College COVID-19 Response Team: Report 9 |
| Śmiertelność | 1.25% | Przyjęte na podstawie Imperial College COVID-19 Response Team: Report 9 |
| Bazowy współczynnik reprodukcji \(R_0\) | 2.0-2.6 | Przyjęte na podstawie Imperial College COVID-19 Response Team: Report 9 |
Tabela 1. Parametry modelu epidemiologicznego COVID-19. Model jest dopasowany do danych przy użyciu podejścia opartego na maksymalnym prawdopodobieństwie, poprzez porównanie przewidywanej liczby zgonów z zaobserwowaną liczbą zgonów. Zakłada się, że dzienna liczba zgonów jest zgodna z rozkładem Poissona. Używając wartości parametrów założonych powyżej, szacujemy wartości \(T_\text{seed}\) oraz \(\kappa\), które najlepiej wyjaśniają zaobserwowane dane przy założonym modelu. Rozważamy cztery różne wartości \(R_0=\{2.0, 2.2, 2.4, 2.6\}\).
Wyniki
0. Kluczowy parametr \(\kappa\) (kappa)
Z punktu widzenia predykcji modelu, parametrem mającym największy wpływ na impakt pandemii COVID-19 w Polsce jest \(\kappa\). Jego wartość odzwierciedla skuteczność interwencji dystansowania społecznego wprowadzonych rzez rząd polski poprzez wpływ na bazowy współczynnik reprodukcji \(R_0\) (reprezentującym ilość drugorzędnych zakażeń przez zakaźną osobę w pełni podatnej populacji; epidemia się rozwija gdy \(R_0>1\) i gaśnie gdy \(R_0<1\)). Konkretniej mówiąc, \(\kappa=1\) odpowiada sytuacji, w której interwencje nie miały żadnego wpływu na transmisję wirusa w populacji (\(R_0\) niezmienione), natomiast \(\kappa=0\) odpowiada sytuacji, w której interwencje uniemożliwiły jakiekolwek nowe zakażenia (\(R_0=0\)). Te obie sytuacje są hipotetyczne i mało prawdopodobne; w rzeczywistości \(\kappa\) będzie przyjmować wartość gdzieś między 1 a 0. Celem zaproponowanego podejścia jest oszacowanie wartości parametru \(\kappa\), jendak jak na dzień 26/03/2020 nie mamy dostatecznej ilości danych, żeby to zrobić. Dla każdej z założonych wartości \(R_0\), szacowana wartość \(\kappa\) była blisko 1, ale wynik ten nie był istotny statystycznie (95% przedział ufności: [0,1]).
1. Kiedy polska epidemia się zaczęła naprawdę?
Z drugiej strony szacunki parametru \(T_\text{seed}\) były istotne statystycznie, co pomogo mi oszacować czas rozpoczęcia epidemii w Polsce (tzn. kiedy epidemia została “zasiana”). Wyniki, pokazane na Ryc. 2, sugerują, że polska epidemia rozpoczęła się gdzieś w drugiej połowie stycznia. Warto odnotować, że to nie znaczy, iż wtedy wirus fizycznie znajdował się na terenach Polski, bo jest wysoce prawdopodobne, że prawdziwy “pacjent zero” zaraził się wirusem za granicą i dopiero po jakimś czasie przywiózł ze sobą infekcję do kraju. Wyniki te natomiast sugerują, że wirus krążył po Polsce na kilka tygodni przed tym, gdy odnotowany został nasz pierwszy przypadek zakażenia koronawirusem 04/03.
Rycina 2. Szacowany czas zasiania epidemii w Polsce. Widełki pokazują 95% przedział nieufności w szacowanym czasie, \(T_\text{seed}\), przetłumaczonych na daty dla zakładanej wartości \(R_0\). Przerywana linia pionowa reprezentuje początek lutego.
2. Ile mamy obecnie zakażonych osób w Polsce?
Jednym z wielkich wyzwań walki z pandemią COVID-19 jest to, że stosunkowo wiele osób przechodzi infekcję bezobjawowo lub łagodnie - przynajmniej przez pewien okres czasu - co utrudnia wykrycie wszystkich zakażonych osób. Niedawny raport z London School of Hygiene and Tropical Medicine sugeruje, że (przynajmniej na dzień 23.03) w niektórych krajach nawet mniej niż 10% wszystkich osób zakażonych SARS-CoV-2 jest oficjalnie zgłaszanych i rejestrowanych. Ilość codziennie wykonywanych testów RT-PCR przeprowadzanych w Polsce zwiększa się z dnia na dzień, ale nie jest jasne, czy coraz lepiej radzimy sobie z takimi niedoszacowaniami, czy też nie. Korzystając z zaproponowanego tu podejścia analitycznego starałem się ocenić skalę i trend w niedoszacowaniu liczby zakażonych, zakładając, że model potrafi dobrze przewidzieć rzeczywistą liczbę przypadków. Wyniki przedstawiono na Rycinie 3. Pokazują one, że w Polsce mamy do czynienia ze znacznym stopniem niedoszacowania liczby zakażonych, jednak jego zmiana w czasie mocno zależy od skuteczności interwencji wprowadzonych przez rząd. Szacuję, że w dniu wprowadzenia tych interwencji, 14/03, zgłoszone przypadki stanowiły około 6%-8% wszystkich zakażonych w Polsce. Dane sugerują jednak, że Polska coraz lepiej radzi sobie z wykrywalnością zakażonych wirusem SARS-CoV-2, ponieważ odsetek ten wzrasta w czasie. Od 26/03, stwierdzam, że wykryliśmy gdzieś pomiędzy 15% do 57% wszystkich zakażonych osób, w zależności od założeń co do wartości \(\kappa\). Szacunek ten nie jest bardzo zaskakujący biorąc pod uwagę epidemiologię choroby i możliwości testowania w wielu krajach. Na przykład wspomniany powyżej raport LSHTM sugeruje, że na dzień 23 marca odsetek wykrywalności liczby zakażonych wahał się między krajami od poniżej 20% (np. Włochy, Hiszpania, Turcja, Wielka Brytania, Francja) do około 50% do 95% (Niemcy, Korea Południowa). Takie porównania należy jednak przeprowadzać z dużą ostrożnością, ponieważ różne kraje znajdują się w różnych fazach epidemii, a skala niedoszacowania powinna z czasem ulec zmianie (mając nadzieję, że się zmniejszy).
Jeśli założymy, że w rzeczywistości wartość \(\kappa\) znajduje się gdzieś pomiędzy 0.2 a 0.8, to w Polsce powinniśmy obecnie wykrywać gdzieś w okolicach 18%-45% wszystkich zakażonych. Co ważne zauważyłem, że trend wzrostowy w liczbie testów wykonywanych codziennie (29700 testów wykonywanych w dniu 26/03) może być niewystarczająco szybki, jeśli transmisja nie została znacząco ograniczona, bo wtedy zakażonych przybywa szybciej niż testowanych. Jak pokazano na Rycinie 3 (panel środkowy), dla wartości \(\kappa\) zbliżonych do 1, stwierdziłem, że liczba testów na zakażoną osobę spada lub maleje, w zależności od wartości \(R_0\). Natomiast dla wartości \(\kappa\) bliższych 0 taka liczba stale rośnie. Biorąc pod uwagę, że \(\kappa\) znajduje się gdzieś pomiędzy 0 a 1, jest prawdopodobne, że wykrywamy coraz większy odsetek wszystkich przypadków w Polsce. Na przykład, jeśli założymy, że \(R_0=2.6\) i \(\kappa=0.5\) (bazowy współczynnik reprodukcji został zmniejszony o połowę do 1.3 przez dystansowanie społeczne), to na 04/03 wykryliśmy (0.4% wszystkich przypadków, na 14/03 wykryliśmy (7.7% wszystkich przypadków, na 26/03 wykryliśmy (28.8% wszystkich przypadków. Niemniej jednak, oszacowanie efektywności interwencji rządowych jest niezbędne do dokładniejszego oszacowania odsetka wykrywalności wszystkich zakażonych w kraju oraz jego zmiany w czasie.