電子回路シミュレータで解く感染症モデル(その2)
〜 実効再生産数Rt を求めよう〜

ABOUTこの記事をかいた人

門 信一郎

門 信一郎

京都大学エネルギー理工学研究所 准教授
佐賀県立伊万里高等学校出身。京都大学理学部卒。九州大学大学院総合理工学研究科修了。博士(工学)。
自然科学研究機構核融合科学研究所助手(現助教)、東京大学高温プラズマ研究センター、東京大学大学院工学系研究科原子力国際専攻准教授を経て2013年2月より現職(宇治キャンパス)。
専門:プラズマ理工学、核融合学、プラズマ計測、分光学。科学教育。
趣味:ピアノ。囲碁。元フィギュアスケート選手。
5回転ジャンプと核融合発電、人類はどちらを先に手にするでしょうか。世代を超えた継続的かつ効率的な育成システム構築が重要でしょう。
神田 峻彦

神田 峻彦

日本工業大学工学部電気電子工学科卒。日本工業大学大学院工学研究科電気工学専攻博士前期課程修了、修士(工学)。電源装置メーカーに勤務、電子回路の設計を担当、現在に至る。
趣味:電子工作、プログラミング。 片道1時間半の通勤時間の有効活用で毎日ノートPCを持ち歩いて、帰りの電車の中で回路図書いたりシミュレーションしたりプログラミングしたり楽しくやってます。月に1回は秋葉原に足を運んで部品を買って、ものづくりに明け暮れてます。休日は月2回、小中学生に科学や電子工作を教えるボランティアをして科学の面白さを伝えています。

目次

1.はじめに
2.感染症モデルの概略 II 「実効再生産数」
2-1) 実効再生産数Rtの「発生間隔計数法」による導出
2-2) 発生間隔計数法の数学的意味と電子回路
3.伝達関数の導入
3-1) 伝達関数による実効再生産数の導出
3-2) 電子回路シミュレータによる伝達関数回路の構築
3-3) 電子回路シミュレータによる実効再生産数の計算結果
4.結語 〜 結果の考え方について
あとがき
参考文献・サイト
付録B:発生間隔、倍加時間、時定数などの関係
付録C:ガンマ関数の「必殺技」

(その1目次へ)


1.はじめに

 

本稿「その1」[1]ではアナログコンピュータ(電子回路シミュレータ)をもちいて、SIRモデルの微分方程式を解くことができることを示しました。主要な素子としては、倍率を与える係数器、2つ以上の電圧の加算値を出力する加算器、入力電圧の時間積分を出力する積分器、2つ以上の入力電圧を乗算する乗算器でした。

本稿「その2」では、電子回路やシステム制御の分野で広く使われる「伝達関数」の考え方に基づいて、実効再生産数を回路出力として表現することに成功いたしました。日々の感染者数(確定日、または発症日)の計数値(人数)を入力しますので、電子回路シミュレータLTspice[2]を入手すれば、どなたでも実効再生産数を求めることができます。ここでの主要な素子は、電子回路部品として、より馴染みのある、抵抗やコンデンサ、オペアンプ(演算増幅器)のみです。電気回路に詳しい方は、もしかすると、より適切な回路構成を構築することができるかもしれません。

「その1」同様、本稿でもちいたソースファイル、および必要な入力ファイル例を公開しています(公開時点)。是非活用いただければ幸いです。

 

2.感染症モデルの概略 II「実効再生産数」

 

2-1) 実効再生産数Rtの「発生間隔計数法」による導出

現在、専門家会議の提言で用いられる評価法も公開されており[3]、前著[4]の分類で言えば、8章の高度なモデリングと統計処理に支えられた「発生間隔計数法」(発病間隔とも訳されているようです。)に基づいた説明が見られるようになってきました。そこでは、単純なSIRモデルにおける回復率(または除去率)を表すγに相当するものとして「発生間隔(Serial Interval)」が現れます。

あらためて発生間隔計数法[4]における実効再生産数の表記を書くと

(接触機会×感染確率) × 感染性を有する時間 = ×1/γ

と表せます。感染性を有する時間は単純なSIRモデル[1]では、[I]にいる時間ですが、「発生間隔計数法」では、感染数の履歴を見ることで、モデルが従うべき「発生間隔」を推定し、その特徴に従って現時点の実効再生産数を推定します。

例えばRt = 1のとき、単純モデルでは「感染数が倍に増える間に、現有感染者が半分になる」と解釈できました。それに対して、発生間隔計数法では「発生数が一定であるとき、2次感染を発生する時定数は、感染性を有する平均時間に一致する」と解釈できます。すなわち[I]は感染性を有する人の仮想的な人口とも言えます(半減期法で、1/γを仮想的な回復の時定数と解釈したのと対照的)。

もう少し数学的に厳密に解釈してみましょう。単純なSIRモデルでは、[I]から[R]への回復・隔離が指数関数exp(-γt)に従うとするものでした。その回復・隔離が指数関数ではなく、ある分布をもった関数g(t)だったらどうでしょうか。そのg(t)に従い、[I]が[I]でなくなります。つまり、感染性を喪失するわけです。すると、[I]と[S]が接触して[S]から[I]に移るレートが変化するので、[I]自身の変化に反映されてくることになります。

いわゆる「西浦モデル」(一般的な統計手法の適用であるが、日本のコロナ禍における北海道大学・西浦博氏の貢献で有名になった。)における感染動向は、上記のような考え方から、以下の「再生産方程式」で表されます[3]。

ここで、J(t)は感染時刻(カレンダー時刻)における新規発病者の計数値であり、g(τ)ははは感染性時間(発生間隔)の確率密度、すなわち、単純なSIRモデルの指数関数的減衰exp(-γt)にかわるものです。τは時間に関する積分変数で、遅れ時間を意味します。

 

感染者数の公表データと合わせるには、「報告日」あるいは「確定日」ではなく「発症日」に補正したデータを用います。これを、「発症日(onset)ベース」の実効再生産数と呼び、明示することにしましょう。もちろん、確定日ごとの計数値を用いれば、確定当日の再生産数として概算値ながらリアルタイムで求めることができます。これを「確定日(confirmed)ベース」の実効再生産数と呼びましょう。ただし、確定日のデータには、無症状感染者も潜伏期間中の検査陽性も含まれるため、解釈には注意が必要です。発症日からさかのぼり、潜伏期間を統計的に補正し、「感染日(infection)ベース」として評価することが求められる場合もあります。

いずれの場合でも、この発生間隔計数法では、いったん適切な発生間隔を推定できれば(ここは高度な疫学統計が必要)、以降、その発生間隔の確率密度を固定し、実効再生産数Rtを推定していきます。施策や治療、隔離によってその関数の形状に多少の変化があったとしても、その変化をRtの値が引き受ける、とも言えるでしょう。

 

再生産方程式(式(1))の具体的なイメージをつかむため、図1に表計算を用いた実効再生産数の導出例を示します。まず、実効再生産数が1、すなわち、感染者が毎日10名ずつ出て、増えない状態にいるとします(図1(a)参照)。二次感染は2日後から1-2-3-2-1-1人と現れるとします。この例では、発生間隔の期待値を計算すると4.3日になります。このとき、たとえば10日後には、過去4〜9日目、つまりJ(t-τ)の遅れ時間 τ として 1〜6日の分 )の一次感染者の割合を足し合わせていることになります。このときg(τ)は全部足して1になるような割合(確率密度)です。過去から見たらその後の発生(表の横列右向き)、現在から見たら過去にさかのぼった影響(表の縦列上向き)ですので、0.1-0.2-0.3-0.2-0.1-0.1となります。

 

図1.発生間隔の期待値(平均値)が4.3日となるような、ある統計分布を仮定した場合の、実効再生産数導出の数値的イメージ(仮定した統計分布に依存する)。実際の感染者を、過去のデータから予測される2次感染者発生数でわったものが実効再生産数Rtに相当する。(a) Rt = 1のとき。(b)日々の増加割合が1.1倍だったとき。

 

もし、感染者の報告が日を追って増えているなら、その割合をかけ合わせたものが、足し合わされることになります。その増え方が実効再生産数Rt に相当します。

そこで、日々の増加率が10%ずつ新規感染者が増加する場合を考えてみます(図1(b)参照)。報告日毎に同じ統計分布g(τ)をかけ合わせてそれを図1(b)の縦に足し合わせて予測値をたてます。実際の新規感染者をその予測値で割ると、実効再生産数1.36が得られます。

単純モデルと対応させると日々1.1倍は、倍加時間7.27日(時定数10.5日)ですので、1/γを回復(隔離)の時定数とすると、半減期(に対応する数値)は1.36×7.27 = 9.9日に相当します。累計の感染者の倍加時間は、割合が一定であれば、数が増えてくれば日々の報告数の増加率と等しくなってきます(割合が変化する場合は注意が必要です)。

一方、発生間隔の期待値4.30を1/γと仮定すると、前著[4]の式(18)より、Rt = 4.3/10.5 + 1 = 1.41 となります。1.36との差はg(τ)の分布の形状に依存するものと解釈できるでしょう。

発生間隔あるいは世代時間と実時間における半減期や倍加時間などの時定数との対応関係を、単純な条件想定のもと、付録Bに考察しました。

 

2-2) 発生間隔計数法の数学的意味と電子回路

このように、この再生産方程式(式(1))は、過去時刻における感染が、次の感染を生む確率を分布g(τ)に従い、現時刻の感染者に「畳み込まれる」ことを意味します。「その人」が感染させるのではなく、統計的な増加の動向を現在から過去を振り返って評価するものともいえます。これを数学的には「convolution(畳み込み積分)」といいます。

ここで、電子回路や線形システム、実験物理学を学んだ方は「装置関数」「インパルス応答」「伝達関数」という用語を耳にしたことがあると思われます。

パルスの入力があっても、検出器で得られる信号は、装置特有の「時定数」によって「なまった」信号になるのは比較的知られた事実でしょう。

電気回路では、この「なまり」を意図的に作ることができます。「ローパスフィルタ(low-pass filter; LPF)」、「積分回路」がそれにあたります。コンデンサによる電荷蓄積やコイルの自己誘導をつかって、パルス時間を遅らせます。大学教養レベルの数学・物理学や電気数学に出てくる「ラプラス変換」をご存知であれば、インパルス応答(時間tの関数)と伝達関数(複素角周波数sの関数)はラプラス変換と逆変換の関係にある、といえばおわかりでしょうか。「伝達関数を求める」とは「インパルス応答を求める」こと、「伝達関数を通す」とは「インパルス応答を畳み込んだ波形を得る」ことと同義と考えてください。

 

図2.インパルス応答(伝達関数)を表す、RCフィルタ回路。

 

例えば、図2(a)の回路では、電源がつくる 5V 0.1秒の矩形の入力パルス(図2(b)、赤線の矩形)に対し、右側にある出力端子ではゆっくり立ち上がり、ゆっくり減衰します(図2(b)、青線のカーブ)。直感的にいうと、連続的な電圧も矩形パルスの組み合わせだと思えば、出力もこのような変形した出力を組み合わせたものになるはずです。これが畳み込み積分の物理的意味です。この単純なRC回路ではこの変化の「時定数τ」は抵抗値R、静電容量Cの積で、

と表されます。図2(a)の回路定数RCで計算するとτ = 100 msで、実際、出力電圧は0.1秒のオーダー(大きさの程度)の変化をしているのがみてとれるでしょう。

 

3.伝達関数の導入

 

3-1)  伝達関数による実効再生産数の導出

本章では、伝達関数を用いて、実効再生産数Rtを求めてみましょう。伝達関数(インパルス応答 g(t))の形状は論文[5]に近いものを作りましたので、厳密ではないにしても、さほど大きく異るものではないでしょう(論文[5]の図1の縦軸は1日あたりではなく、0.1日あたりの数値のようです)。 採用したg(t)具体的な形は、大学での数学や統計学でよく知られるガンマ関数Γ(*)を用いた以下の式(3)になります。

ただし、n = 1、m = 3、1/a = 1.1  です(aを逆数で書くのは、時定数を表すため考えやすいからです)。この形を採用する利点や計算方法は付録Cに述べています。

いく通りかの関数型を図3に示します。

 

図3.発生間隔(serial interval)の確率密度。n, m, aはそれぞれ、次数、重み、時定数の逆数、に関連した式(3)における係数。これらの組み合わせで、適切な関数型を構成することを試みた。

式(1)の右辺の時間積分は伝達関数(インパルス応答g(τ))が0になると効かなくなりますので、2週間程度でも十分でしょう(τ無限大まで積分しても結果は同じ)。この部分は伝達関数の回路を通すことで自動的にやってくれます。

あとは、日々の発症者報告数と、回路の出力値との比をとります。図4にその概念図を示します。

 

図4.電子回路をもちいて発症者ベースで実効再生産数Rtの求める概念図。平滑化は目的と必要に応じて行う。

 

伝達関数を求めるには、図1の横の列と同じ考え方で、1人を1日間だけ入力したときに、何日後にどのような割合で感染を広げるかを示す確率分布を作ることを目標にします。その確率分布はすでに評価済みの(論文等に出て認められた)もの(例えば文献 [5])、あるいはそれに類する形状のもの(例えば式(3))にします。そのような伝達関数が得られると、日々の発症者数の入力に対して、その人たちが次の感染を起こす人数の予想値が得られます(図1の縦列)。

データのばらつきが大きい場合、データをそのまま用いて演算処理(特に割り算や微分)を行うと、そのばらつきがますます増大し、データの解釈ができなくなることがよくあります。したがって、実験データを扱う場合は適切な「前処理」をして次のステップに進むことが行われています。平均操作もその一つです。時間的に平均をすることと、何度も実験をやって、回数分の平均をとること(アンサンブル平均)は、通常、同等とみなされます。

データを滑らかにする代表的な方法は「平滑化(スムージング)」と言われるもので、例えば、ある1点のデータの前後数点を平均し、それをその1点の代表値として表現します。次の点は、またその前後の平均値をとります。濃い霧の中のスポットライトの移動をイメージしてみてください。1点の光が滑らかに動くかわりに、その光はその周囲の場所にも影響を与えることになります。どの程度の「幅」で平均をとるか、あるいは、その「幅」のなかに強弱をつけるかどうかは、データの品質と解析の目的とのバランスを保つように選びます。一連の同様な解析には、同じ処理の仕方で比較検討を行うことが重要です。

本稿では、平滑化回路として、図2に示したような「ローパスフィルタ」を用いています。電気回路で、高周波成分をカットするために広く使われており、ノイズキャンセラーの一方式として有名かもしれません。どの程度の周波数(時間変化の速さ)をカットするか、というのは、これも、データの品質と目的に応じて決定します。

 

3-2)  電子回路シミュレータによる伝達関数回路の構築

では、これをLTspiceで実装してみましょう。回路図を図5(a)に示します。

図5.電子回路シミュレータを使った、実効再生産数導出のための(a)伝達関数回路 (b)伝達関数への入力(単位インパルス波形)および出力波形(インパルス応答g(t))。

伝達関数の回路は、非反転増幅の2次ローパスフィルタとして代表的なものです。各RCの値を探索するには、「その1」でおこなったSIRモデル構築と同様、粒子群最適化を用いました。伝達関数(インパルス応答)は、時間積分が入力と出力で同じになるように規格化されたものです(図5(b))。入力は「デルタ関数」[6]を実効的に模擬するため、100V-0.01秒のパルスにしました(単位インパルスと言います)。規格化は出力前のR5、R6の抵抗比を用いて調整します。今の値は1.01でしたので、誤差1%です。これは得られる精度に対し、無視できる量とみなせます。

オペアンプなどの電子回路に不慣れな方は、図2のような単純なRC回路を多段に接続して、近い形状の伝達関数を作成してみるところからチャレンジしてみられてはどうでしょうか。「感触」がつかめるのではないかと思います(実用上はなかなか適度な時間に0に落とすのがむずかしいですが、ピーク位置をあわせ、時間積分を1に近づけることを目標に)。

本回路で採用した伝達関数(インパルス応答g(t))は、実際に専門家会議の評価で採用されている再生産方程式(式(1))における発生間隔の確率密度g(τ)や、それを模した数式(3)のg(t)と比べると若干の違いはあります。これを合わせ込むには、フィルター回路の構造を高次にするなどで改善が見込めますが、本稿ではあえて典型的な回路のまま近づける試みをしました。実際、用いるデータのばらつきや誤差の影響もあり、決定される実効再生産数は、伝達関数の形状の変化には鈍感といえます。

 

3-3)  電子回路シミュレータによる実効再生産数の計算結果

 

図6.電子回路シミュレータを使ってもとめた日本における発症日ベースの実効再生産数の推移。入力は発症時間がわかっているもののみの計数を用いた(17593例中、9026例)。6月以降は発症日がまだわからないデータがあることを考慮し表示範囲から省いてある。入力の起点は2020年1月3日であるが、「.tran」コマンドで、記録開始を30日後に設定したため、表示上の0秒=2月2日である。

 

計算結果を図6に示します。実効再生産数Rtは、入力電圧(実際の新規発症者数)を出力電圧(Rt = 1 の場合の発症者の予測値)で割ったものです。これもLTspice内で計算して表示し数値データとして書き出すことができます。発症日毎の計数データは、「新型コロナウイルス感染者数マップ」[7]にまとめられている感染者リストのファイルからカウントしました。発症日データのない報告もあるのですが、本稿では原理実証のため、入手可能な発症日のデータをそのまま用いました。入力する計数として、誤差を許容し、データ数と速報性を重視すれば、確定日のものを選ぶと、確定日ベースの実効再生産数が得られます。都道府県別の計数値を用いればその県の実効再生産数が得られます。

結果は、専門家会議等[8][9]から公表されている実効再生産数の時間的推移をほぼ反映しており、評価値もかなり近い値が得られているのがわかります。発症者ピークの伝達関数への入力波形と出力波形の中間は4月4日付近です。ちょうどそのあたりで実効再生産数はRt = 1を横切っています。ただし、政府の発表にあるような「感染日ベース」の実効再生産数に比べると、「発症日ベース」のものはピークが5日程度遅れ、「確定日ベース」のものは、さらに5日程度遅れます。(次章図8参照)。

濃厚接触者を追跡し、通知するアプリケーション等が浸透し、データが蓄積してくれば、もしかすると感染日ベースの実効再生産数も直接評価が可能になるかもしれません。

 

実効再生産数が少し「バタついて」いるのは、入力データのバタつきを反映しています。入力する前に数値的に移動平均などの前処理(平滑化・スムージング等)をすることによって、なめらかになります。ただしその分、時間分解能は低下しますので、もとのデータの性質にあわせて適切な方法をとるのが通例です。図6に示した発症日ベースの実効再生産数の評価では、平滑化をせずに伝達関数に入力しています。「確定日」を用いると検査体制のばらつきを反映し、さらにバタつきますので、フィルタ処理が必要だと思われます。

 

平滑化も回路シミュレータを用いる場合、図4のように、適当な時定数の平滑化回路を前段に挿入することで入力信号をなめらかにすることが可能です。図2に示したRC回路も代表的な平滑化回路の一つです。

実際に、東京都の確定日ごとの計数値を用いて、平滑化回路を通し、確定日ベースの実効再生産数を求めてみましょう。結果を図7に示します。この場合は、データのばらつきが大きいので、前段に式(2)で与えられる時定数として4秒(すなわちグラフでは4日)の平滑化回路を入れています。

 

図7.(a)平滑化回路を組み込んだ伝達関数の回路図。(b) 電子回路シミュレータを使って求めた東京都のにおける確定日ベースの実効再生産数の推移。東京は発症日データがほとんど得られていない(2020年6月29日データで確定日6172中、発症日の記載は258人のみ) 。入力の起点は2020年1月3日であるが、「.tran」コマンドで、記録開始を60日後に設定したため、表示上の0秒=3月3日である。時間軸を7日刻みに表示すると、曜日毎の傾向がわかりやすい。

 

東洋経済オンライン[10]で公開されている実効再生産数とほぼ同程度時間スケールで相対的な挙動を示しています。このサイトで採用されている、西浦氏の提案による簡易式(付録Bの式(B5))では1週間分の総計を求めているため、当日のデータが、おおよそ1週間の中央、すなわち3.5日前の移動平均値となっています。一方、図7(b)については、生データと平滑化後の出力をみると3日前後の遅れがあるのがわかります。その原因は平滑化回路の時定数(式(2))によるものです。生データを、数値的に平滑化し、それを伝達関数に入力すれば、この遅れをなくすことができるので、より早い段階で新規感染者計数値の変化を検知することが可能です。

平滑化回路へ入力すると、この信号遅延のため、出力値が最新日よりも先まで得られます。そこで、平滑化回路の入力と出力の傾向が重なる程度まで、伝達関数の入力(=平滑化回路出力)と出力の時間軸を手前にずらせば、回路シミュレータのみで数値的平滑化と同様の結果が得られます。本稿では回路出力の特徴を示すため、あえてそのまま表示しています。

このように、回路シミュレータを利用すると、データ処理のどの段階で、データのばらつきや遅れが評価に響いてくるかを確認しながらデータ解析を進めることができる、という利点もあります。

 

4.結語 〜 結果の考え方について

 

現在にいたるまでの経時の計数データをもとに、現在の(実際には、発症数が確定した日における)実効再生産数を求めることは、感染が拡大しているのか、収束しているのか、その目安を把握するために意義があると思います。種々の評価方法が公開されており、発生間隔の関数型と日々の計数値からRtを求めるエクセルシートが適用されている論文も出ています[11][12]。ただし、その数学的な「中身」を理解するにはやや敷居が高いようにも思えます。

回路シミュレータを用いれば、いったん回路を構築してしまえば、あるいは動作が確認された回路を利用できれば、数式や統計処理をブラックボックスとして、実効再生産数の概算値、その経時変化を追っていくことができることがわかりました。こちらも一つの手法として提案できる方法ではないかと思います。複数の手法で同一の結果が得られることは、「SIRモデル」や「実効再生産数」などの概念の信頼性向上や様々な分野の人たちへの浸透にも有用と思われます。

もちろん、もし厳密性を要求するならば、個々人について、図8に示すように、

 

感染潜伏発症  → 受診・検査確定・報告  治癒・死亡 

 

といった、一連の時間的な状態の変化を加味する必要があります。それぞれ個人差があるので、統計分布(ポアソン過程など)に従う、などの統計的処理を行います。

 

図8.(a)感染から発症、確定診断を経て治癒または死亡に至るまでの時間のばらつきの模式図。本稿で用いたデータは17594名中、9026名の発症日の計数を用いた。発症から確定までかかった日数のヒストグラムを(b)に示す。

 

ですが、施策の決定が細かい数値の差に依存するような場になければ、信頼区間の評価など、数値そのものに付帯する幅にある程度目をつぶり、迅速性と、データの入手性を考慮して相対変化による傾向を把握することも一つの考え方だと思います。

この個人差、特に感染から発症までの誤差は小さいと仮定することにより、その日の発症者の数がほぼ確定した段階で、リアルタイムに実効再生産数を出すことができます。もっとも、発症から確定までには1週間程度の遅れはあると思いますし、感染は発症の数日前です。やはり行動変容が、数字として現れるには2週間程度の遅れがあるのは致し方ありません。

かといって、隔離のために無症状の人たちに対する大規模な検査を定期的に実施するのは、クラスター追跡の過程や興行団体など公共性が高い集団の場合を除き、現実的ではないと思います。将来的には、濃厚接触の機会を検知するスマートフォン用のアプリケーションを用い、通知された人への漏れのないフォローアップ検査が徹底できれば、感染拡大に対する予測性の向上にも有効かもしれません。ですが、その有効性が現れるには、体制整備などに期間を要すると思われます。

通常は行動変容を維持しつつ、症状がある場合は自主的に隔離、マスクなど人に移さない防護措置をすること、発症者について迅速に検査、判定を行うことによって、発症日と確定日のタイムラグを最小限に抑えることで、感染の拡大を抑制し、実効再生産数の評価の遅れも改善されると思われます。

 

今般のコロナ禍は本稿執筆の2020年6月26日現在、いまだ進行中であり、日本よりも深刻な事態に陥ってしまった国も多くあります。日本も全く気を抜けない状況です。今後の「災害復興」で必要になることは、恐らく日々発表される感染状況のデータに一喜一憂することなく、「社会的・経済的要因を考慮に入れながら合理的に達成できる限り」リスクを下げながら、感染拡大防止に努めつつ、数字やその変化の動向を冷静に判断できるようになることではないでしょうか。

 

科学的なリテラシーを基礎として一人ひとりが主体的に現状を理解し、施策の意味や効果を冷静に判断できるような社会の醸成、そのために本稿が少しでもお役に立つことを、著者としても心より願っております。

 


 

あとがき 

 

門:本稿の執筆の機会をくださいました京都女子大学の水野義之名誉教授(@y_mizuno)に感謝いたします。SNSで、LTspiceを用いてSIRモデルを解かれている[twitter_1]の投稿に感銘をうけて、(お名前も存じ上げなかった)神田氏(@Akibajin)に突然お声をかけて議論を進めました。同氏には、私の様々な提案(無茶振り)にも適切に対応いただき、内容がどんどん深まっていきました。得意分野の異なる者同士が気軽に連絡をとれるのもSNSの魅力でもあると実感しました。

 

神田:「感染症の流行をコンピュータ上で模擬することができないか」と興味本位で調べていくとSIRモデルにたどり着きました。数理モデルが微分方程式で表されていたので、アナログ回路で演算できないかと思い、回路シミュレータで試行錯誤したのが始まりでした。普段からSNSで趣味の電子工作関連の発信をしているので、いつもの流れで発信したところ[twitter_2]や[twitter_3]など、思わぬ反響とともに門先生、京都女子大学の水野先生との出会いがありました。

無茶振り(?)はありましたが、さまざまなアドバイスをいただき、とても充実した結果を得ることができました。今回このような共著という形で貴重な機会をいただき、心から感謝するとともに、厚く御礼申し上げます。

 


 

参考文献・サイト:

 

[1] 門 信一郎、神田 峻彦:電子回路シミュレータで解く感染症モデル(その1)〜SIRモデルを解こう〜 , RAD-IT21 WEBマガジン (2020.07.14)  https://rad-it21.com/サイエンス/kado_kanda_20200714_1/

[2] アナログ・デバイセズ社。LTspiceダウンロードサイト。https://www.analog.com/jp/design-center/design-tools-and-calculators/ltspice-simulator.html#

[3]西浦博 GitHubリポジトリ   (2020.5.12閲覧) https://github.com/contactmodel/COVID19-Japan-Reff

「JASTJ COVID-19 科学ジャーナリストのための情報整理」で補足解説。 https://note.com/jastj/n/n463f1aeb1833

[4] 門 信一郎:この感染は拡大か収束か:再生産数Rの物理的意味と決定, RAD-IT21 WEBマガジン (2020.03.27)  https://rad-it21.com/サイエンス/kado-shinichiro_20200327/

[5] H. Nishiura et al., “Serial interval of novel coronavirus (COVID-19) infections” International Journal of Infectious Diseases 93 (2020) 284–286. https://doi.org/10.1016/j.ijid.2020.02.060

[6] デルタ関数: https://ja.wikipedia.org/wiki/ディラックのデルタ関数 (2020年6月22日閲覧)

[7] 「新型コロナウイルス感染者数マップ」について、(2020.6.20時点) https://jag-japan.com/covid19map-readme/

[8]「新型コロナウイルス感染症対策の状況分析・提言」(2020年4月22日)https://www.mhlw.go.jp/stf/seisakunitsuite/bunya/0000121431_00093.html

[9] 新型コロナウイルス感染症対策専門家会議: https://ja.wikipedia.org/wiki/新型コロナウイルス感染症対策専門家会議   (2020年6月22日閲覧)

[10] 東洋経済オンライン特設ページ(新型コロナウイルス国内感染の状況)https://toyokeizai.net/sp/visual/tko/covid19/

[11] 山中伸弥による新型コロナウイルス情報発信「実効再生産数(Rt)の計算を試みました」https://www.covid19-yamanaka.com/index.html (6/20閲覧) 内に紹介された文献; Anne Cori 他, American Journal of Epidemiology,178,(2013)1505-1512. https://doi.org/10.1093/aje/kwt133 に発生間隔の関数型と日々の計数値からRtを求めるエクセルシートがある。(2020年6月22日閲覧)

[12] Rt Covid-19 Japan, https://rt-live-japan.com, およびその引用 (2020年6月22日閲覧)

 


 

付録B 発生間隔、倍加時間、時定数などの関係

 

実効再生産数Rtはしばしば「感染症の流行が進行中の集団のある時刻における、1人の感染者が生み出した二次感染者数の平均値」と説明されます。文字通りに解釈するとRtの累乗で新規の感染者が増え続けることになります。

ではこれを実時間軸にあわせ、単純モデルと対応させることができるでしょうか。仮に、ある時間間隔τがすぎたときに一斉に感染が起きる特殊な場合を想定して計算してみましょう。

最初1名だったとすると、現有感染者数I(t)は

と表せます。新規感染者が出た時点で、それまでの感染者は感染性をなくすからです。

一方、SIRモデルでは、感染初期において、倍加時間をτ2とすると、

なので、両者より

と関連づけられそうです。この特殊な想定下でRt =2であれば、τは倍加時間τ2と一致します(実際は統計分布に従いますので「なまり」ます)。その意味で、「発生間隔毎の等比級数的な感染拡大」は、「現有感染者の半減期の間に倍加時間が何回入るぐらいの感染拡大かを表す」、という単純モデル(半減期法)によるRtの物理的意味と合致します。よくつかうRtの値(0.5〜6、ただし1を除く)であれば、両時間の比τ/τ2 は概ね0.5倍から2倍程度の値をとります。

この関係を使うと、西浦氏が[3]にて簡易方式として提案し、東洋経済のサイト[10]で評価されている実効再生産数Rtの速報値の意味が納得できます。左辺はτ日毎にRt倍、右辺は7日毎に新規感染者がX倍と書くと、現有感染者(新規感染者数に等しい)の時間変化は式(B1)より

となります。発生間隔の平均値は約5日なので、τ = 5を代入し、

が得られます。この簡易式を確定日(報告日)で整理したデータに用いると、過去2週間の確定人数を使って当日の実効再生産数が求まります(発症日のデータに適用できれば発症日ベースの実効再生産数が得られる)。よって、得られる値は半週程度前の実効再生産数概算値を反映しているといえます。さらに、7日分の計数値にはすでにそれ以前の伝達関数が畳み込まれているので、その分の誤差も含みます。したがって、実際のRtの急な変化に対し、評価値はやや「鈍感」にならざるを得ません。

本稿の伝達関数を利用する方法では、同じ確定数の経時データを用いると、伝達関数(式(3)または図3、図5(b))の広がりは、およそ9-10日程度です。したがって、クラスターの発生や、行動変容が実効再生産数Rtの値に現れるまでの時間を、4−5日程度短縮することができると考えられます。
 
 

参考:K値について

同様の考え方で2020年4月中旬から5月にかけて発表[B1]されて以来、話題になっているK値とRtとの関係を評価することができます。K値の定義は

で与えられています。ここで1週間毎に新規感染者がX倍に増えるなら、初項1からn週過ぎた時、分子はXn-1、の累計感染者数は初項1の等比級数の和で

です。ただし、統計分布を考慮しない簡易式(B4)よりX= Rt (7/5)です。

この方程式はおそらく解析的には解けませんので、解く場合は数値的に行います。ですが、式(B7)からわかるように、統計分布を考慮しないRtの簡易値である式(B5)とK値は対応しています(Rtが一定であれば、Rtが決まればK値が決まる)。K値は0〜1の値をとり、収束トレンドで0に近い値になります。どちらも経時変化の傾向の指標になりますが、筆者(門)は、K値が感染者累計数、すなわち感染規模に依存するので、収束トレンドの直感的なイメージはつかみにくいのではないかと思っています。例えば、日々の感染者が一定数である状況が続いている場合を考えてみます。実効再生産数はRt = 1で一定ですが、K値は分子が1のまま、分母が等差数列の和n(n+1)/2に従い増え続けますので、nが増えるにしたがい、0に向かって減少を続けます。そこからわずかにRtが増えれば感染爆発ですが、K値は減り続けます。K値の「減り方の変化」を検知する必要がありますが、それはRtの変化をみるのと同じです。しかも、K値の場合はどこまで値を戻せばよいかが得られません(分母が異なっているから)。したがって、どこかでカウンタをリセットして計数を再開することになります。K値の「リセット」は、トレンドの可視化と感度の向上に有用という言い方もできます。しかしその反面、値には任意性が避けられないことも分かります。

 

[B1] 中野貴志教授(大阪大学核物理研究センター)による論文等(K値について)

 https://www.osaka-u.ac.jp/ja/news/info/corona/corona_info/from_members/rcnp_nakano

 


 

■付録C: ガンマ関数の「必殺技」

 

伝達関数(インパルス応答g(t))は積分すると1になるように「規格化」しないといけませんが、式(3)の積分にはガンマ関数Γ(n)[C1]というものを導入した、便利な公式が使えます。mnを整数、aを正の実数とすると式(C1)が得られ、式(C2)の公式を用いて簡単にすることができます。

ただし、n!は 整数nの階乗を表します。

大学で学ぶ物理学において、量子力学に現れる波動関数の計算ではn = 1を、気体分子運動論、統計力学に現れるマックスウェル・ボルツマン分布を用いた計算ではn = 2をよく使います。証明はガンマ関数の定義から変数変換して簡単に求められます。筆者(門)が大学生の頃、この種の計算(ガウス積分とも言われ、複素積分や漸化式で導出する、、、と習う)に嫌気がさして変数変換して以来、「必殺技」と呼んで重宝しています。

式(3)は「必殺技」をもちいると、つねに積分が1になるようにできていることがわかります。

[C1] https://ja.wikipedia.org/wiki/ガンマ関数   (2020.6.20閲覧)

[C2] https://ja.wikipedia.org/wiki/ガンマ分布 (2020.6.20閲覧)

 


 

(2020年6月22日投稿、7月14日公開)

 
 

ABOUTこの記事をかいた人

門 信一郎

門 信一郎

京都大学エネルギー理工学研究所 准教授
佐賀県立伊万里高等学校出身。京都大学理学部卒。九州大学大学院総合理工学研究科修了。博士(工学)。
自然科学研究機構核融合科学研究所助手(現助教)、東京大学高温プラズマ研究センター、東京大学大学院工学系研究科原子力国際専攻准教授を経て2013年2月より現職(宇治キャンパス)。
専門:プラズマ理工学、核融合学、プラズマ計測、分光学。科学教育。
趣味:ピアノ。囲碁。元フィギュアスケート選手。
5回転ジャンプと核融合発電、人類はどちらを先に手にするでしょうか。世代を超えた継続的かつ効率的な育成システム構築が重要でしょう。
神田 峻彦

神田 峻彦

日本工業大学工学部電気電子工学科卒。日本工業大学大学院工学研究科電気工学専攻博士前期課程修了、修士(工学)。電源装置メーカーに勤務、電子回路の設計を担当、現在に至る。
趣味:電子工作、プログラミング。 片道1時間半の通勤時間の有効活用で毎日ノートPCを持ち歩いて、帰りの電車の中で回路図書いたりシミュレーションしたりプログラミングしたり楽しくやってます。月に1回は秋葉原に足を運んで部品を買って、ものづくりに明け暮れてます。休日は月2回、小中学生に科学や電子工作を教えるボランティアをして科学の面白さを伝えています。