グラフ関数電卓で観る感染症数理モデル

2020/05/06
ダウンロードリンク修正 2020/05/08
SIR2.g1m 差替えアップデート 2020/05/09
プログラムアップデート 2020/05/11

グラフ関数電卓で、感染症数理モデル - SIRモデル - をグラフ化

新型コロナウィルスのパンデミックに際して、底の浅い情報、信憑性が疑われる情報、噂が蔓延しています。これをインフォデミックと言うそうです。新型コロナウィルスについては、まだよく分かっておらず未知のウィルスと言うべき段階なので、自分の分かる範囲で情報の取捨選択をして納得できる理解をし、それを逐次修正して行動に反映したいと思っています。

地上波TVの情報はあてにならないと感じており、主にインターネットやYouTubeで情報を集めているなかで、感染対策の考え方に東大方式と京大方式があったり、大阪モデルのように自治体によって明確に異なるものがあったり、台湾方式やスウェーデン方式が対局にあるのですが国ごとに違いがあることを知りました。

様々な違いがある一方で、世界各国の対策全てに共通する感染症の数理モデルの基本が、およそ１００年前のスペイン風邪流行期に提出された SIRモデルだと知りました。この数理モデルを少しでも理解すれば、実際にパンデミック対策の基本の考え方を知る良い機会になりそうです。そこで、先ずSIRモデルをグラフにしで、次に重要なパラメータを変えたらグラフがどのように変化するかをビジュアル化しようと思います。

SIRモデルの微分方程式は難しいので、差分形式に変更して数値的にアプローチすれば、エクセルでグラフ化できると考え、その具体的な計算方法を調べました。そして、エクセルでの計算方法をグラフ関数電卓の Casio Basic に移植しました。

※実際のグラフ描画のコードは Colon様の支援を頂きました。グラフ描画について経験が殆ど無い管理人には、とても良い勉強になりました。

  
fx-CGシリーズ　　　　　　　　　　　　　　　　　fx-9860シリーズ

高精細度カラー液晶のfx-CGシリーズの方がグラフが圧倒的に観やすいですね。

グラフを描く時に使ったパラメータは、以下の "パラメータ一覧・入力画面" で、変更してグラフ化できます。
 

これらパラメータの意味については、さらに読み進んでもらえば、お分かり頂けると思います。

⇒ ファイルのダウンロード
ダウンロードした SIR_Model2.zip を展開すると、純正Casio Basic用プログラムファイルは、g3mフォルダとg1mフォルダにあります。

▶ fx-CGシリーズ：g3mフォルダにある下記4つのファイルを全て fx-CGシリーズの PROGRAMフォルダに転送し、純正Casio Basicで SIR2 を実行してください。
　・SIR2.g3m
　・INPI.g3m（入力ボックス 2.1G - 整数版）
　・INP.g3m（入力ボックス 2.1G - 小数版）
　・3DS.g3m（3桁区切り）

▶ fx-9860Gシリーズ  (fx-9860G OS2.00 以降のモデル)：g1mフォルダにある下記4つのファイルを全て fx-9860GシリーズのPROGRAMフォルダに転送し、純正Casio Basicで SIR2 を実行してください。
　・SIR2.g1m
　・INPI.g1m（入力ボックス 2.1G - 整数版）
　・INP.g1m（入力ボックス 2.1G - 分数版）
　・3DS.g1m（3桁区切り）
※ fx-9860G OS1.xx では MENUコマンドが実装されていないため、OS2.00 以降で実行可能
※ fx-9860G OS2.00 以降には、fx-9860GII、fx-9860GIII を含む
※ このプログラムは、現行の C.Baisc for CG Ver 1.45 や C.Basic for FX Ver 2.54 では、純正Casio Basicのグラフ描画プログラムの全てにまだ対応していないので、動作しません。


[2020/05/11: ver 2.1 に差替えアップデート]

• グラフデータを格納したリストの最大値を求める関数を変更してスピードアップしました
• グラフデータ準備に多少の時間がかかるので、その間 = Processing data = と表示するように変更しました

▋目　次
１．はじめに
２．ＳＩＲモデルと再生産数 - 接触機会の抑制 80% の謎
３．ＳＩＲモデルの可視化 - エクセルでグラフを描く
４．純正 Casio Basic への移植
５．おわりに

---

1. はじめに

新型コロナウィルス (COVID-19) の世界的な感染拡大で、私たちの生活が大きく変わってしまっています。ご自身やご家族の健康的被害、そして経済的被害に遭われている方々には心からお見舞い申し上げます。さらに、医療関係者、そして国民の生活に必須な生産や物流に携わっている方々に、深い敬意と感謝を捧げたいと思います。

日本は、感染症に対する医療的リソース (専門家と施設) が他の先進諸国に比べて圧倒的に不足している状況で、また強制力を伴った措置の法的立て付けが無いにも関わらず、死者数が欧米諸国に比べて1桁から2桁以上少ない、といったかなり特殊な状況にあります。

2020年5月6日の時点では全体的に感染者数拡大が鈍化する傾向があり、医療施設や介護施設、そして家庭内の感染が顕著になっている状況です。緊急自体宣言が5月31日まで延長されましたが、未知の新型コロナウィルスを "正しく恐れる" ことで、個人が自分の行動を工夫することで、早期に収束が迎えられるのではないか、是非そうなって欲しい、と思います。

施設内感染や家庭内感染の調査により、物に付着しているウィルスを人が触れて感染する傾向が高いことが分かってきています。このモノーヒト感染においては、ウィルスが自分の口や鼻、目に入らないことが重要で、ウィルスが付着しやすい手で顔を触らず頻繁に消毒することがかなり有効な感染対策とされています。

マスクは、自分から他人、そして他人から自分への飛沫感染を抑制するだけでなく、口と鼻を手で触らないためにも有効です。手で顔を触らない習慣は身につくものだそうです。さらに頻繁に石鹸での手洗いやアルコールでの消毒を行えば、モノーヒト感染はかなり効果的に抑えられます。

ヒトーヒト感染については、人と人の物理的接触によるウィルスの移動、そして飛沫感染を抑える効果を狙っていると考えられます。日本政府は "接触機会の抑制" を 80% にすることを要請しています。１人の外出を 55% 抑制すると、２人が出会う確率は 0.45 × 0.45 = 0.2025 ≒ 20% になり、３人以上ではさらに低い数値になります。外出を 55% 抑制すれば 80% の接触機会の抑制が達成されます。

新型コロナウィルスは、空中のウィルスがマスクをすり抜けて感染するような "空気感染" は殆どないらしいことも分かってきています。

従って、マスクをして、手で顔を触らず、手洗いや消毒を実施し、外出を 55% 抑制すれば、感染収束に向かう可能性が高いのではないでしょうか？

今は未知のウィルスを相手にしているのですから、ワクチンも治療薬も無く、万一発症すれば医療機関は命を救うためには、出来うる限りの対処療法を施すしかありません。そのための医療関係スタッフと医療機器、そして感染症用の施設、つまり感染症対策の医療リソースを、余裕を持って確保することが最重要課題となっています。

そのためには、感染防止対策を行いながら、感染者数の動向をリアルタイムで追いかけ、今後の動向予測が不可欠です。未知の新型コロナウィルスなので、統計的、疫学的に動向をシミュレーションして、医療リソースが不足しないように施策をうつことが最優先です。そこでは、感染症の数理モデルのシミュレーションが拠り所になるのでしょう。

インフルエンザでは患者数の動向を追いかけるのに、新型コロナウィルスでは何故感染者数で追いかけているのか、不思議に思っています。未知のウィルスだから数理モデルを用いて "感染者" で整理するのは判りますが、患者数、重傷者数、死者数の動向も重要なはずだと思っています。

感染症や医療のプロでなくとも、我々一般国民は論理的、科学的思考を常にアップデートしてゆきたいものです。そして未知の感染症に対して "正しく恐れる" 態度を持ちたいと思っています。

▲目次に戻る

---

2．ＳＩＲモデルと再生産数 - 接触機会の抑制 80％ の謎

ＳＩＲモデルは、ほぼ100年前の1920～1930年の間にマルケックとマケンドリックにより提出された数理モデルです。このモデルの基本的な考え方は、今日に至るまで使われているので、それだけ信用に足るものと言えます。ネットで調べると、SIRモデルを基点として様々な発展型モデルが提案され、実際に利用されていることも分かります。

 
これが、SIRモデルの常微分方程式です。

式を見ているだけでは、なんだかよく分からないのですが、このモデルから導き出される重要な指標の1つが、1人の感染者が何人に感染させるかの指標 -  "再生産数" です。再生産数には、R₀、R_e、R_t があります。

• R₀：基本再生産数。何も対策をとらない時の病原体の素の感染力を示す。
• R_e：有効再生産数。R₀に手洗い、マスク、行動制限などの対策後の感染力を示す。
• R_t：実行再生産数。感染拡大防止対策に加えて集団免疫率を加味した時の感染力を示す。

基本再生産数 R₀ は、英語では R nought (Rノート) や R Zero と呼ばれます。2011年公開の映画 "コンテイジョン" では R nought といっています（最近観なおしました）。R₀ は数理モデルの分析から推定される値で、新型コロナウィルスではまだ確定的な値が分かっていません。現在のところWHOによれば、新型コロナウィルスの R₀ は、1.4～2.5 程度の可能性が考えられているようです。

▶ 代表的感染症のR₀：(国立感染症研究所の資料による）

• はしか：12-18（飛沫感染）
• 風疹：5-7（飛沫感染）
• おたふく風邪：4-7（飛沫感染）
• ポリオ：5-7（経口感染）
• 天然痘：6-7（飛沫、接触感染）
• インフルエンザ(スペイン風邪)：2-3（飛沫感染）
• ジフテリア：6-7（飛沫感染）
• 百日咳：12-17（空気感染）

実効再生産数 R_t は、感染拡大対策の要因に加えて集団免疫率を加味したものだと思われます。免疫を獲得する人は感染症が拡大するにつれ増えるので、実行再生産数 R_t は感染症の拡大と収束に伴って変化してくる値です。ドイツなどでは R_t を公表して、国民への指針の根拠としています。合理的ですね。日本でも最近になって 実行再生産数 R_t という言葉を使うようになっています。

さて、日本政府が使っている "接触機会の抑制" という言葉は、どうも不明瞭、不明確に感じています。
ここでは勝手に、
"接触機会の抑制率" = "ウィルスに接触する機会の抑制" = "接触抑制率"
と定義します。

ここで定義した "接触抑制率"を p としたとき、この指標だけで考える有効再生産数は、R_e = R₀ × (1 - p) となります。

R_e は累乗で効いてくるので、R_e = R₀ × (1 - p) < 1 の時は、感染拡大の抑制効果が高くなります。但し、０でない限り感染拡大初期では感染が広がります。これはＳＩＲモデルの式から分かります。

新型コロナウィルスの基本再生産数 R₀ は、専門家でもまだよく分かっていないので、例えば大げさに大きな値として 2.5 ～ 4.8 までの値になる時、R_e < 1 とするための接触抑制率は、この式から分かります。 

 
多分 R₀ が 4.8 ということは無い (もっと小さい) でしょうが、あり得ないくらい悪い状況だとしても、接触機会を 0.8 つまり 80% 抑制すれば、R_e が1未満になるわけです。

新型コロナウィルスの R₀ は 1.4 ～ 2.5 程度だろうとWHO (国際保健機関) が言っているので、この範囲での最悪の値 2.5 の場合は、接触抑制率を 80％ にすれば R_e は 0.5 となり、拡大防止策に多少手抜かりがあっても余裕を持って感染を抑えられそうですね。おそらくこのあたりが、現在の指針である 80％ 抑制の背景ではないかと思われます。

▲目次に戻る

---

3．ＳＩＲモデルの可視化 - エクセルでグラフを描く

常微分方程式を素直に解いて、式を導くのは難しいですが、グラフを描くのが目的なので、差分形式に変換すればなんとかなると考えました。そこで、下記のSIRモデルの常微分方程式から、
 

と、少し工夫した上で、差分形式に変形してみました。



ここで、パラメータとして、

• R₀：基本再生産数
• p：接触抑制率 (上述)
• b：接触抑制を行った時の実際の再生産数 (有効再生産数、Re)
• c：一旦感染した人が感染力を失う（免疫を獲得するか死亡する）割合
• N：感染が始まる前 (t=0) の人口 = S(0) の値
• I(Δt)：感染が始まる最初 (t=Δt) の感染者数
• Δt：差分の刻み

とすると、時刻 t+Δt における S(t+Δt)、I(t+Δt)、R(t+Δt) の値を S(t)、I(t)、R(t) を使った漸化式で計算できそうです。

そこで、エクセルには以下のように計算式を入れました（クリックすると拡大できます）。


計算式は150行もあれば十分でしょうか？

散布図を描いてみると、以下のようになります。

実際には、右端にある黒い背景のセル（5箇所）は、表で対応するセルへリンクしています。

⇒ ファイルのダウンロード
ダウンロードした SIR_Model2.zip にエクセルファイルを同梱していので、これをご覧ください。

パラメータを変えると、赤い感染者数のピークの高さや位置が変化する様子が分かります。
基本再生産数 R_o や、接触抑制率 p を変えて 有効再生産数 R_e が大きく1を下回ると、感染者のピークが見えなくなるくらいに小さくなり、感染が拡大しないことがグラフで可視化できます。

最初の感染人口を多くすると、ピークの位置が左にシフトして感染拡大が早まることが分かりますが、感染者ピークの高さはさほど変わりません。これはチョット意外です。

現実に感染者が増え、医療機関のお世話になる人が増えすぎると、医療機関がパンクします。感染者数が医療機関の能力を超えないようにしたいので、有効再生産数 R_e をどの程度コントロールすれば良いのかがポイントになります。

R_e を小さくすると、感染者ピークの高さが小さくなり、ピークの幅が広がり、ピークの位置が右にずれます。これが対策の狙いだということが、グラフの変化からよく理解できます。

※ 微分液式から差分形式への変換は、100％の自信がありませんので、批判的にご覧頂き、間違いがあればご指摘下さい。

▲目次に戻る

---

4．純正Casio Basicへの移植

さて、上の差分方程式が正しいものとして、これをCasio Basicに移植します。

実際のグラフ描画の部分は、当ブログの強力サポータである Colon様の手によるものです。これが SIR (Ver 1) です。
管理人は、差分化形式での計算方法を Colon様に提供し、パラメータ一覧・入力部分の作成を行い、SIR2 (Ver 2) にしました。

管理人は、グラフ描画関係の関数をこれまでプログラムで使ったことが無かったので、Colon様のコードはとても勉強になります。きちんと動くプログラムのソースはをよく読むのは、最善の勉強方法ですね。SIR2 のグラフ描画部分は殆ど SIR と同じです。


▍プログラムの基本構成

t、S(t)、R(t)、I(t) の値をエクセルでの計算式に基づいて、それぞれ List 1、List 2、List 3、List 4 に格納し、このリストで 数式を定義し、最後に DrawFTG-Con 関数か DrawFTG-Plt 関数を用いて、定義した数式をグラフとして描画します。
なお、SIR2 では時間刻みΔt を 1 に固定しました。

fx-CG50 (OS 3.30) / fx-CG20 (OS 3.11) のソフトウェアマニュアル、8-31 ページに以下の記載があります；


これだけでは簡潔すぎて分かりにくいと思いますので、プログラムの詳細については実際に SIR2 のソースを見て下さい。


▍プログラムを使ってみる

SIR2 を起動すると、パラメータ一覧・入力画面になります。以下のように初期値を決めました。
　・[0]キーでグラフ描画
　・[1]キーで基本再生産数 R₀ を変更
　・[2]キーで接触抑制率 p を変更
　・[3]キーで隔離率(回復率) c を変更
　・[4]キーで最初の人口 N を変更
　・[5]キーで時刻1での初期感染者数を変更
　・[6]キーで横軸(時間軸)の最後を変更
　・[EXIT]キーでプログラム終了
　※ 有効再生産数 R_e も表示しています

新型コロナウィルスに罹患してから回復するのに２週間程度かかることが分かってきましたので、１週間あたりの回復率＝隔離率は 0.5 として良さそうです。これを初期値にしました。それに関連して、時間刻みを１週間としました。

パラメータ一覧・入力画面で、[0]キーを押すとメニューが現れます。
 

ここで、Connect タイプか Plot タイプを選択して [EXE] を押すとグラフが描画されます。

  
Connectタイプ　　　　　　　　　　　　　　　　  Plotタイプ

感染者のピークは 28週に現れます。

さて、このグラフ描画で使った初期設定のパラメータでは、有効再生産数 R_e = 1 となっています。これは１人が1人を感染させる状況ですので、感染は広がります。

この状態で [F1]を押すと、3つのグラフのどれかを [↑] / [↓] キーで選択できます。赤い色の式が表示されたら、赤い感染者数のグラフが選択されたことになり、分かりやすいです。そして、[←] / [→] キーでグラフに沿ってカーソルを移動でき、カーソルのある位置の座標値が表示されます。これで、ピークの位置にカーソルを移動すると、そこが28周目だと分かります。

これは、グラフ関数電卓内蔵のグラフ描画コマンドに内蔵されている機能です。

グラフが表示された状態で [EXE]を押すとメニューが表示されるので、そこで 3:Reset Parameters を選びます。
 
ここで、[EXE] か [3] を押すとパラメータ一覧・入力画面に戻るので、
接触抑制率 p を 0.7、つまり70%に変更してから、[0]を押してグラフを再び描画します。

   
有効再生産数 R_e が 0.75 と 1 未満になり、赤い感染者ピークが半分以下になり、ピークの位置が47週まで遅れることが分かります。医療リソースに余裕が生まれ、必要な体制を整えるための時間が稼げます。


再びパラメータ一覧･入力画面に戻ります。ここで 接触抑制率 p  をさらに上げて 0.75 にして、グラフを描画します。

  
有効再生産数 R_e ＝0.62 とさらに低くなると、感染者ピークが60週まででは見えにくくなっています。

そこで、ｘ軸（時間軸）を 90 週まで長くするように変更してみます。

  
すると、感染者ピークが大幅に減り、 72週まで大幅に遅れることが分かります。
接触抑制率を上げるのは、確かに大きな効果があることがグラフで可視化できました。

有効再生産数 R_e を下げる効果が確認できました。
実行再生産数 R_t は集団免疫率までも加味したものなので、R_e よりは小さい値になると思います。

最も単純なSIRモデルに、接触抑制率 p だけを含めた今回の極めて単純なモデルでも、おおまかな傾向が分かりました。


▍メモリエラーへの対処
ところで、List に多くのデータを格納するので、残りメモリが少ないとメモリ不足によるエラーが発生することがあります。
 

この場合は、メインのアイコンメニューで Memory を選び、メモリマネージャーでリストファイルを削除すれば問題は解決します。
  
メモリマネージャで[F1]を押して、F1:Main Memory を選び、

  

LISTFILES フォルダを [F1]で SELECT し、[F6]でDELETE します。すると削除確認のポップアップが現れるので、Yes:[F1] を押します。

  

すると、リストファイルがフォルダごと削除されたことが確認できます。これでメモリの問題は解決です。

▲目次に戻る

---

５．おわりに

最初の１年は R_e が小さくなるように行動制限をかけ、次の１年は行動制限を緩くして R_e が大きくなったとき、感染者ピークのグラフがどのように変化するかを調べるシミュレーションにするなど、SIR2 プログラムには色々な改造案が考えられます。

今後の政府の方針を見ながら、現実に即したシミュレーションプログラムを作ってみては如何でしょうか？

皆様からのプログラムのご投稿やご提案を歓迎致します。

▲目次に戻る

---

▶ 管理人の独り言：

未知の感染症への対応は、疫学的な正当性と経済的な正当性のせめぎ合いだと思います。疫学的正当性を貫けば経済的死者数（自殺）が増え、経済的正当性を貫けば病死者数が増えます。このせめぎ合いのなかで、死者数を最小化するための着地点を見いださなければなりません。

現在の未知の新型コロナウィルスの蔓延は、過去に事例がありません。特に最近の日本では、幸いなことに感染症(SARS、MARS、鳥インフルエンザなど) の被害があまり大きくならずに済んでいて、諸外国に比べて感染症に関する日本国民全体の危機感が低かったと思います。すなわち国民の代表である政治家も危機感を持っていなかったと思います。

過去に例をみない事態を処理するのは、霞ヶ関の官僚や地方のお役人には苦手な分野です。平時には法律に則って緻密な行政を執り行うのが役人の本来業務です。法律を超えることは許されていません。前例に捕らわれず、新たな法律を作り、時には超法規的な措置を行うのは政治家の役割、政治家にしかできないことです。政治家と役人が緊密に連携し、現実を科学的に理解し、大胆な発想で国民の生命や財産を守って頂きたいと思います。

明治時代の先達は、西欧の Economics に "経済" という日本語を与えてくれました。この言葉は「世を治め、民を救う」という意味の "経世済民" が語源です。疫学的な解決だけでなく、バランスの良い疫学対策と経済対策を強力に進めて頂くよう切に願います。

我々国民は、政治家のパフォーマンスをしっかりチェックし、我々の権利である投票で正しく判断したいものです。


当ブログでは異例の意見表明ですが、今回だけは例外としてお許しください。





応援クリックをお願いします。励みになるので．．．

にほんブログ村

人気ブログランキング

FC2ブログランキングへ

 


keywords: プログラム関数電卓、コロナウィルス グラフ関数電卓

---

リンク集　|　ブログ内マップ