Pandemi COVID-19 melalui mata ahli matematika, atau mengapa model SEIRD klasik tidak berfungsi

Abstrak, atau tentang waktu luang para ilmuwan muda


Selama beberapa minggu terakhir, rekan saya dan saya telah mengakhiri hari kerja dengan bersaing dalam akurasi perkiraan untuk pengembangan epidemi COVID-19 di Rusia menggunakan berbagai metode regresi non-linear. Dan jika ramalan untuk hari esok ternyata menjadi baik, maka ramalan untuk periode lebih dari satu minggu hanya mencerminkan kenyataan secara umum. Tampaknya semuanya sudah jelas: ada model epidemiologis , ada metode optimasi , ada data yang cukup rinci , cukup untuk menggabungkan ini bersama-sama dan mendapatkan perkiraan yang akurat selama sebulan, atau bahkan enam bulan, di muka. Pada artikel ini saya akan membagikan pemikiran saya tentang apa yang salah dengan model SEIRD klasik dan bagaimana cara memperbaikinya. Dan, tentu saja, saya akan membuka tabir kerahasiaan yang membungkus masa depan kita dengan Anda.

Duduk, seorang iblis yang marah menunggu kita untuk mereka yang tahu apa persamaan diferensial itu (untuk yang lain, gambar-gambar indah terlampir).


Gambar di atas menunjukkan jumlah total kasus COVID-19 yang dikonfirmasi pada skala logaritmik untuk Rusia dan tiga negara Eropa dalam 5 besar dalam hal jumlah yang terinfeksi. Penjelasan lebih lanjut dalam teks.


Menit Perawatan UFO


Pandemi COVID-19, infeksi pernafasan akut yang berpotensi parah yang disebabkan oleh coronavirus SARS-CoV-2 (2019-nCoV), telah secara resmi diumumkan di dunia. Ada banyak informasi tentang Habré tentang topik ini - selalu ingat bahwa Habré dapat diandalkan / bermanfaat, dan sebaliknya.

Kami mendesak Anda untuk kritis terhadap informasi apa pun yang dipublikasikan.


Sumber resmi

, .

Cuci tangan, rawat orang yang Anda cintai, tinggal di rumah kapan saja memungkinkan dan bekerja dari jarak jauh.

Baca publikasi tentang: coronavirus | kerja jarak jauh

Model SEIRD


Model epidemi SEIRD milik kelas yang disebut model kompartemen , yang intinya adalah untuk membagi populasi menjadi beberapa kelompok ( kompartemen bahasa Inggris ), dalam kasus kami:$ S $( Bahasa Inggris rentan) - rentan,$ E $( Eng. Terkena) - mereka yang memiliki penyakit pada masa inkubasi,$ I $( Bahasa Inggris menular) - sakit,$ R $( Bahasa Inggris pulih) - pulih,$ D $( Inggris mati) - orang mati. Kemudian, ukuran masing-masing kelompok dibandingkan dengan variabel dalam sistem persamaan diferensial, penyelesaiannya, Anda dapat memprediksi dinamika epidemi. Ada banyak modifikasi dari model SEIRD, misalnya, SEIR adalah model yang disederhanakan yang tidak secara terpisah mempertimbangkan kasus-kasus pemulihan dan kematian. Untuk berkenalan dengan model lain, saya bisa merekomendasikan artikel yang bagus tentang topik ini.

Sedikit teori


Untuk pertama kalinya, model epidemi dalam bentuk sistem tiga persamaan diferensial untuk variabel $ S, I, R $muncul dalam karya W. Kermak dan A. McKendrick pada tahun 1927.
Persamaan diferensial ini memiliki bentuk:

$ \ begin {align} \ frac {dS} {dt} & = - \ beta \ frac {SI} {N}, \\ \ frac {dI} {dt} & = \ beta \ frac {SI} {N} - \ gamma I, \\ \ frac {dR} {dt} & = \ gamma I, \ end {align} $


di mana, di samping variabel yang kita kenal, konstanta berikut muncul: $ N = S + I + R $ - ukuran total populasi, $ \ beta $ - tingkat penularan infeksi, $ \ gamma $- kecepatan pemulihan.

Arti persamaan Kermak dan McKendrick adalah sebagai berikut: jumlah kerentanan menurun secara proporsional dengan jumlah mereka dikalikan proporsi rata-rata yang terinfeksi dalam populasi.$ I / n $, jumlah yang terinfeksi tumbuh pada kecepatan yang sama, disesuaikan dengan kenyataan bahwa beberapa dari mereka $ \ gamma I $pulih, dan jumlah pemulihan meningkat karena penurunan jumlah yang terinfeksi. Perlu dicatat bahwa model SIR mengandung nonlinier$ SI $, karena itu solusi analitik dari sistem persamaan menjadi umumnya tidak mungkin, tetapi, untungnya, metode diferensiasi numerik dapat dengan mudah mengatasi tugas ini.

Menambahkan variabel lain di sini$ E $ (jumlah orang dengan penyakit pada masa inkubasi), kami mendapatkan model SEIR:

$ \ begin {align} \ 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} $


di mana konstanta lain muncul $ \ kappa $- laju transisi penyakit dari tahap inkubasi ke kondisi terbuka. Gambar dari diambil dari artikel .



Segera jelas bahwa model SEIR tidak sangat cocok untuk menggambarkan COVID-19, jika hanya karena dalam model ini ada pembawa infeksi yang tersembunyi.$ E $tidak menular. Kekurangan ini dapat diperbaiki dengan memperkenalkan, mengikuti Pengpeng et al. parameter opsional$ \ theta $, mengkarakterisasi tingkat infektivitas pembawa infeksi laten dibandingkan dengan yang sakit. Model SEIR yang dimodifikasi, yang akan kami coba terapkan pada epidemi saat ini, akan terlihat seperti:

$ \ 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} $


Pada pandangan pertama, model yang dihasilkan berjanji akan sepenuhnya dapat dipercaya.

Eksperimen numerik dengan model SEIR


Untuk pemodelan, kami akan mencoba untuk mengambil parameter berikut, dengan fokus pada data terbuka . Dengan asumsi bahwa penyakit berlangsung rata-rata 14 hari (setidaknya berapa lama bentuk ringan berlangsung , yang menyumbang hingga 80% dari kasus), kami menemukan nilainya$ \ gamma = 1/14 = $ 0,0714. Akan menerima$ \ beta = 3/14 = $ 0,2143. Besarnya$ \ theta = $ 0,6dipinjam dari Pengpeng et al . Mengingat masa inkubasi rata-rata 3 hari, ambil$ \ kappa = 1/3 = $ 0.33. Kami menerima populasi Rusia yang sama$ N = 144,5 \ cdot10 ^ 6 $orang.

Sebagai kondisi awal, kami menggunakan data untuk Rusia pada 2 April, ketika langkah-langkah yang diperkenalkan pada akhir Maret untuk membatasi penyebaran infeksi memiliki efeknya, yaitu:

$\begin{align} &S_0=3548,\\ &I_0= 3283,\\ &E_0=0,5 I_0. \end{align}$


Peringkat $E_0$kami mengambilnya secara relatif sewenang-wenang, dan itu tidak masalah, karena, seperti yang Anda tahu, ada sesuatu yang salah.

Sebagai hasil pemodelan dengan metode Euler dengan langkah 1 hari dari 2 April hingga 24 April inklusif, kami memperoleh grafik yang ditunjukkan di bawah ini: di sebelah kiri dalam skala linier, di sebelah kanan dalam logaritmik.



Penanda bulat menandai data nyata tentang jumlah total kasus di Rusia, penanda persegi pada jumlah pasien. Sekilas, hasilnya terlihat bagus, kecuali satu: dengan parameter model, kami jelas tidak menebak. Dan di sini metode optimasi datang membantu kami.

Optimalkan


Metode optimasi, jika pembaca tidak terbiasa dengan mereka, adalah algoritma yang memungkinkan menemukan minimum beberapa fungsi obyektif. Dalam kasus kami, kami miliki sebelum kita masalah regresi nonlinier: bagaimana memilih vektor parameter persamaan diferensial$\mathbf{x} = (\beta, \gamma, \kappa, E_0)^\top$ sehingga himpunan titik solusi dari persamaan diferensial $F$ sedekat mungkin dengan set poin pengamatan $\mathcal{F}$.

Kami menggunakan standar deviasi sebagai ukuran kesalahan model. Fungsi objektif akan berbentuk

$ 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}, $


Dimana $M$ - jumlah poin $F$ - jumlah total kasus infeksi yang diberikan oleh model, $\mathcal{F}$ - jumlah total kasus aktual, $G$ - jumlah pasien saat ini, yang memberikan model, $\mathcal{G}$- total nyata dari kasus aktif saat ini.

Menggunakan Kotak Alat Pengoptimalan di MATLAB, kami akan menyesuaikan parameter model dengan data pengamatan. Sebagai hasilnya, kami mendapatkan solusi yang ditunjukkan pada gambar di bawah ini.

image

Sekilas, semuanya baik-baik saja. Perbedaan itu ternyata sama$f(\mathbf{x}) = 131,98$, dan "pada mata laut yang cembung", kesesuaian solusinya berjalan dengan baik. Mari kita lihat parameter yang diterima:

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


Nilai hampir 8 juta pasien laten, dengan sekitar 60 ribu kasus tercatat
pada 24 April, adalah sesuatu yang meragukan. Kami juga menemukan bahwa waktu transisi rata-rata ke fase aktif penyakit adalah$1/\kappa = 2079$hari.

Kenapa ini terjadi? Semuanya akan menjadi jelas jika kita menganalisis bentuk kurva pada skala waktu yang lama. Untuk melakukan ini, ambil model SEIR kami dengan parameter "masuk akal" dan simulasikan untuk jangka waktu yang lama (dalam percobaan ini saya mengambil makna baru$\beta=0,186$):



Kurva yang sesuai dengan jumlah kasus memiliki bentuk-S karakteristik pada skala linier. Program optimisasi mencoba memberikan bentuk kurva ini. Selain itu, ramalan itu sendiri dengan parameter "masuk akal" mengerikan - menurutnya, hampir 90% dari populasi negara itu akan sakit pada bulan September - jelas tidak realistis jika Anda melihat hasil untuk negara lain (gambar yang sama seperti pada awal artikel, hanya pada skala linier) ):



Di sini saya membandingkan tiga negara Eropa dalam 5 besar dalam jumlah kasus, dan Rusia. Dapat dilihat bahwa kita berada sekitar satu bulan di belakang laju epidemi, dan bahwa selama sebulan, pertumbuhan jumlah total kasus di ketiga negara hampir linier (dan bahkan lebih lambat dari linier), berbeda dengan hasil yang diperoleh dalam model SEIR. Ini menimbulkan tiga pertanyaan:

  1. ?
  2. SEIR, ?
  3. , , ?

Saya akan mulai dengan menjawab pertanyaan ketiga. Ketika kami memperkirakan sesuatu, kami menghadapi tugas yang agak tidak menyenangkan: data yang kami bangun modelnya tidak ideal - mereka mengandung kesalahan, kebisingan, dan model yang dibangun atas dasar mereka juga mengandung beberapa kesalahan. Ketika kami melanjutkan deret waktu dengan titik-titik model kami, kesalahan terakumulasi - dan agak cepat, jika kami memprediksi fungsi yang meningkat dalam waktu. Dan ini adalah kasus dalam kasus kami. Apalagi model tersebut adalah model yang mencerminkan situasi nyata yang sangat terbatas. Perkembangan epidemi yang mendadak di kota besar baru, penggunaan metode pengobatan yang lebih efektif, perubahan dalam metode pengumpulan informasi - semua ini dapat membuat begitu banyak kesalahan dalam data nyata sehingga perkiraan jangka panjang akan sepenuhnya jauh dari kenyataan.

Kruk dan sepeda: kami memodifikasi model SEIR


Mari kita coba menjawab pertanyaan mengapa pertumbuhan epidemi melambat menjadi linear. Dengan jumlah orang yang terinfeksi yang kita miliki sekarang, peran signifikan dimainkan oleh skala ekonomi yang terkait dengan kecepatan komunikasi yang terbatas antara orang-orang.

Lebih tepatnya, kita ingat: jumlah kasus dalam model SEIR berbanding lurus dengan rata-rata jumlah kasus dalam populasi$I/N$. Aturan ini berfungsi baik dalam populasi kecil, di mana setiap orang dapat berkomunikasi dengan semua orang, dan pasien didistribusikan secara merata. Pada kenyataannya, terutama pada skala puluhan dan ratusan ribu orang, jika Anda mengambil dua orang yang sakit secara acak, ternyata mereka tidak hanya tidak pernah berkomunikasi satu sama lain dan tidak melihat satu sama lain, mereka bahkan tidak naik mobil subway yang sama. Dan memang mereka tinggal di kota yang berbeda. Semua yang menyatukan mereka adalah rantai ikatan sosial, yang mengarah pada fakta bahwa virus itu ditularkan kepada mereka.

Sebagai contoh, saya membuat model epidemi dalam bentuk otomat seluler, di mana setiap sel berinteraksi dengan hanya 4 sel yang berdekatan. Ini setara dengan fakta bahwa setiap individu dalam populasi memiliki 4 kontak sosial - ini adalah jumlah yang sangat kecil untuk populasi manusia, tetapi semakin cepat efek membatasi koneksi sosial akan terwujud dengan sendirinya. Pada setiap iterasi dengan probabilitas 0,1, masing-masing dari 4 tetangga sel yang terinfeksi dapat terinfeksi. Penyakit ini berlangsung rata-rata 14 hari. Hasil simulasi untuk kumpulan 200x200 sel disajikan pada gambar di bawah ini, di mana$k$Apakah nomor iterasi.



Warna biru menunjukkan kerentanan, sakit kuning, hijau pulih. Yang paling menarik adalah seperti apa grafik dari jumlah kasus. Dan mereka tampak seperti yang direncanakan: setelah fase pendek pertumbuhan sub-eksponensial, seperti dalam model SEIR, ada fase pertumbuhan linear yang berlarut-larut - seperti dalam kenyataan.



Tujuan saya bukan untuk mendapatkan gambaran yang menyerupai kenyataan secara kuantitatif. Jika Anda ingin kredibilitas lebih tinggi, saya dapat merekomendasikan proyek Sergey Potekhin, yang baru-baru ini diterbitkan di Habré . Untuk pembaca yang teliti, bukti lebih kuat dari pertumbuhan linear diberikan di bawah ini.

Bukti teorema epidemi pertumbuhan linear dalam populasi besar
: $d$- . . , 20 ( ) $d \approx 4$. $d$- . $n^{\frac{1}{d}}$, , $P$, $2P$. , ,

$n_{k+1}=\left(n_k^{\frac{1}{d}} + 2P \right)^d,$



$n$ : $n_k = n(t_k); n_{k+1} = n(t_k + 1)$. , :

$n'_{k+1}=n'_{k}\left(1 + \frac{2P}{n_k^{\frac{1}{d}}} \right)^{d-1}.$



, $n_k$ $n'_{k+1}=n'_{k}$, , .

4% 16- .



Beralih ke statistik global tentang epidemi, kita akan melihat hal yang sama: selama sebulan sekarang, pertumbuhannya sudah linear, terlepas dari kenyataan bahwa model klasik menjanjikan kelanjutan dari kenaikan eksponensial dalam jumlah kasus.



Tugas untuk yang penasaran
  1. COVID-19, offline .
  2. .1 ?
  3. .1 2, .


Sekarang tentang modifikasi model SEIR. Hal paling sederhana yang dapat kita lakukan adalah memperbanyak komponen non-linear dengan beberapa fungsi, tergantung pada jumlah pasien. Dengan jumlah kasus yang kecil, fungsi ini harus mendekati 1, dengan jumlah yang besar, ia seharusnya cenderung nol. Kandidat yang cocok paling sederhana adalah

$\mathcal{\varphi}(I,E) = e^{-\alpha(I +\theta E)^{K_0}}.$


Pemilihan parameter $\alpha$ dan $K_0$dapat mengimbangi pertumbuhan eksponensial dalam model asli.

Tambahkan ke model, untuk informasi lebih lanjut, dan komponen$D$- jumlah kematian. Dapatkan versi modifikasi dari model SEIRD:

$ \ 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 I. \ end { sejajarkan} $



Hasil simulasi ditunjukkan pada gambar di bawah ini.



Kesalahan standar dibandingkan dengan model asli hampir tidak berubah. Nilai parameter sudah realistis. Untuk kenyamanan, saya menetapkan jumlah awal yang terinfeksi pada fase aktif sebagai$ I_0 $.

$ \ 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} $



Model ini menginterpolasi turunannya dengan sangat baik - nilai peningkatan harian dalam jumlah kasus dan jumlah kematian.



Mari kita coba membuat perkiraan. Kami mengambil cakrawala perkiraan 2 bulan dan terus memodelkan solusi dengan parameter yang ditemukan oleh program optimisasi.



Pada pandangan pertama, tidak buruk, tetapi Anda tidak bisa berharap ramalan seperti itu ke tanah air tercinta Anda: jumlah kasus baru akan terus menurun, tetapi jumlah total akan terus bertambah. Dalam hal ini, epidemi dapat dihentikan hanya dengan bantuan vaksin, atau dengan menunggu sampai hampir seluruh populasi sakit. Jumlah kematian baru ditetapkan sekitar 200 per hari. Ini adalah ilustrasi yang jelas tentang apa yang akan terjadi jika kita tidak memperkuat langkah-langkah untuk memerangi epidemi. Apakah itu yang menanti kita? Dan demi masa depan yang tidak terlalu cerah ini, banyak dari kita yang rajin duduk di rumah, membeli soba dan kertas toilet?

Di bawah ini saya akan mempertimbangkan dua skenario, dan melihat jarak berkabut bulan-bulan mendatang 28 April 2020, saya tidak bisa mengatakan dengan pasti mana yang akan mengembangkan acara lebih lanjut. Sekarang, pada saat memecah kurva kasus-kasus baru, kita berada pada titik di mana ada dua masalah untuk memprediksi sesuatu.

Skrip AS


Hegemon dunia berada dalam posisi yang sangat tidak bisa dihidupkan. Tertunda dalam membuat keputusan kunci yang akan memperlambat pertumbuhan epidemi pada awalnya, ia masih tidak dapat mengatasi peningkatan alami dalam kasus-kasus baru.

Model SEIRD yang dimodifikasi, dilatih pada 33 poin pertama, mulai dari 2 Maret, plus atau minus secara realistis memprediksi perjalanan epidemi pada bulan April.



Seperti yang Anda lihat, pertumbuhan pada bulan April hampir sepenuhnya linier. Model sedikit melebih-lebihkan angka kematian pada bulan April, tetapi gambaran keseluruhannya benar.



Gambar ini menunjukkan peningkatan harian dalam kasus baru dan kematian di Amerika Serikat. Ini sangat mirip dengan apa yang diprediksi model untuk Rusia.

Skenario Jerman


Orang-orang Jerman yang berdisiplin berhasil mengubah kurva demi mereka, dan pertumbuhannya lebih lambat daripada linear. Selain itu, untuk membuat model yang relevan, saya harus secara manual menambahkan peningkatan tingkat pemulihan pada 6 April$ \ gamma $1,7 kali, jika tidak demikian penurunan jumlah kasus yang tajam dalam hal model SEIRD tidak dapat dijelaskan.



Model itu dilatih pada 27 poin pertama, mulai dari 10 Maret. Saya juga mengubah fungsi nonlinear. Untuk Jerman, eksponen yang tergantung waktu lebih baik:

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


Jenis fungsi ini menunjukkan peningkatan kumulatif dalam jumlah ikatan sosial yang terputus dan, dengan demikian, penyebaran infeksi. Di sini Anda memiliki ilustrasi yang jelas tentang manfaat isolasi diri.



Di atas menunjukkan peningkatan harian dalam kasus baru dan kematian. Seperti dalam kasus Amerika Serikat, data nyata mengandung fluktuasi yang diucapkan dengan jangka waktu 7 hari. Ini berarti bahwa pada akhir pekan jumlah kontak meningkat, dan akibatnya, jumlah orang yang terinfeksi bertambah. [UPD: sebaliknya, ini menurun, karena indeks isolasi-diri Yandex dapat secara tidak langsung memberikan kesaksian . ]

Kesimpulan


Membuat prediksi - baik jangka pendek dan jangka panjang - bukan hanya penghargaan untuk rasa ingin tahu. Jika terjadi epidemi, Anda perlu tahu berapa tempat tidur yang harus disiapkan, berapa banyak ventilator yang akan diproduksi, berapa bulan untuk menimbun peralatan pelindung pribadi untuk dokter. Pejabat harus memahami apakah tindakan yang diambil sudah memadai atau apakah larangan dan pembatasan baru harus diberlakukan. Idealnya, model tersebut harus mencerminkan kenyataan dengan sangat baik sehingga memungkinkan untuk melihat kekuatan tindakan dari setiap tindakan yang baru diadopsi, dan kemudian dimungkinkan untuk memperkuat langkah-langkah yang bermanfaat dan membatalkan keputusan yang ternyata tidak berguna.

Meskipun dengan beberapa pemesanan, setiap orang yang belum mendapatkan kekebalan dari COVID-19 dapat benar-benar dijelaskan dengan surat itu$ S $, seluruh variasi orang sakit - surat itu $ I $, dll. Selain itu, model SEIRD bahkan dapat menjelaskan sesuatu. Tapi dia bisa memprediksi sesuatu di masa depan dengan sangat kasar.

Saya sengaja mengutip dalam artikel itu hanya skenario negatif ketika kita mengulangi nasib Amerika Serikat dalam hal dinamika epidemi. Jika skenario ini ternyata benar, pada akhir Juni kita akan memiliki lebih dari 300 ribu kasus penyakit yang terdaftar dan lebih dari 10 ribu kematian. Meskipun ada prasyarat untuk fakta bahwa skenario ini tidak terwujud, saya akan menyarankan Anda untuk menghubungkannya sesuai dengan prinsip: "berharap untuk yang terbaik, bersiaplah untuk yang terburuk." Seperti kata pepatah, jika pikiran terbaik NASA telah berperang melawan epidemi di Amerika Serikat , maka ini benar-benar hal yang buruk.

Sejauh ini, yang tersisa bagi kita adalah mengunjungi tempat-tempat umum lebih sedikit, untuk digunakandengan respirator yang tepat , cuci tangan, jangan lupa untuk menyeka smartphone dengan alkohol setelah Anda mengeluarkannya di jalan, dan ikuti rekomendasi sederhana lainnya.

Tapi tetap saja, mungkinkah membuat prediksi yang lebih akurat? Ya tentu saja. Tetapi tentang ini lain waktu.

Jika Anda memiliki keinginan untuk bermain dengan kode sumber dan menyarankan skenario pengembangan sendiri, maka di sini adalah tautan ke github .

All Articles