What's New

SIRモデル*1

下式が感染症の数理モデルの基本であるSIRモデルの連立微分方程式で、時間(t)の関数はそれぞれ

の各グループを表します。

\( \frac{ \mathrm{ d } S(t) }{ \mathrm{ d } t } =- \beta S(t)I(t) \\ \frac{ \mathrm{ d } I(t) }{ \mathrm{ d } t } = \beta S(t)I(t) - \gamma I(t) \\ \frac{ \mathrm{ d } R(t) }{ \mathrm{ d } t } = \gamma I(t) \)

2番目の式を変形して

\( \frac{ \mathrm{ d } I(t) }{ \mathrm{ d } t } =( \beta S(t) - \gamma )I(t) \)

感染が流行、つまり \( \frac{ \mathrm{ d } I(t) }{ \mathrm{ d } t } \gt 0 \) になるのは

\( \frac{ \beta S(t)}{ \gamma} \gt 1 \)

のときで、また初期条件

\( S(0)=1 \) から

\( \frac{ \beta }{ \gamma} =1\)

が流行の閾値で β γ は基本再生産数(R0)と定義されます。

SIRモデルの連立微分方程式には \(S(t)\) と \(I(t)\) の積の項が含まれます。この積の項がなければ、現代制御理論の状態方程式のように線形な連立微分方程式として解法できそうですが、 \( S(t)I(t) \) の項があるSIRモデルでは、 \( \beta \) と \( \gamma \) のパラメータを適当に設定した上でシミュレーションという手段をとらせざるを得ないようです。

アナログコンピュータ

コンピュータといえばデジタルそのもので、その対極としての古臭さく非効率なアナログというような認識からは、 アナログコンピュータという名称それ自体が論理矛盾のように受け止められるかもしれませんが、今から半世紀前、 パソコンどころかマイコン(マイクロプロセッサや組み込み用のマイクロコントローラ)すら存在していない時代に、アナログコンピュータという装置が確かに存在していました。

実際のアナログコンピュータはどのようなものであったかというと、演算増幅器(オペアンプ)、抵抗(ポテンションメータ)、コンデンサで構成された積分器や加算器等を組み合わせた回路によって微分方程式や制御系の伝達関数をシミュレーションする装置で、工業系の学校で教育用にも使用されていたようですが、本来の目的はリアルタイムな大砲の弾道計算であったという説もあります。

アナログ半導体の代表格オペアンプ(OP-Amp)の名称はOperational Amplifier(演算増幅器)の略称なのですが、元々の由来はアナログコンピュータにあるのではないでしょうか。オペアンプのアプリケーションにステート・バリアブル・フィルタとかバイクワッド・フィルタと呼ばれる回路があります。これらの回路では積分器や加算器を組み合わせて希望する系の伝達関数を実現しており、ある意味極小のアナログコンピュータといえるかもしれません。


SPICEアナログ回路シミュレータ

SPICEとはWindows以前から存在するアナログ回路シミュレータのアプリレーション・ソフトウェアで、MS-DOSではPSpice、UNIXにはHSPICEがありました。Windowsの時代になってもデファクト・スタンダードはPSpiceでしたが、アナログ半導体メーカ旧リニア・テクノロジーからLTspiceが無償提供されるとアマチュアの間でのデファクトになりました。LTspiceは現在でも旧リニア・テクノロジーを買収したアナログ・デバイセズのサイトから無償ダウンロード可能です。

SPICEはデジタル演算でアナログ回路シミュレーションを行いますが、SPICE上で積分器や加算器を組み合わせてアナログコンピュータのように微分方程式をシミュレーションをすることも可能です。 ちょっと倒錯しているかとも思えますが、かなり有効な手法です。 この場合、実際のオペアンプのSPICEモデルは用いず、ゲインを106(120dB)程度に設定した電圧制御電圧源(Voltage Controled Voltage Sorurce)を演算増幅器の代わりに使用します。SPICEモデルとしてのVCVSはゲインを無限大に設定できないことを除いては、入力インピーダンス無限大、出力インピーダンスがゼロ、周波数帯域が無限大の理想オペアンプの条件を満たしています。

このあとアナログ回路シミュレータのLTspice上で感染症SIRモデル微分方程式のシミュレーションに挑戦します。


乗算器をLTspice上で構成

\(S(t)I(t)\) の演算のためにLOGアンプ→加算器→アンチLOGアンプの構成で乗算器を作成しました。抵抗値は実際の電気電子回路としてどうなのか?という値になっていますが、広いダイナミックレンジを確保するためのあくまでシミュレータ上の値ということで割り切ります。

下のグラフは上のシミュレーション回路のDCスイープ出力です。広ダイナミックレンジな二乗特性が得られました。
回路図右上のR6の値はこのグラフをにらみながら入力1V→出力1Vになるようカット＆トライで決定しました。


LTspiceによるSIRモデルシミュレーション

下の回路図はLTspiceによるSIRモデルシミュレーションの全回路です。

下のグラフは上回路でのシミュレーション結果です。V(st),V(it),V(rt) の他に乗算器の動作確認するため、乗算器出力 V(stxst) と StとItそれぞれの出力同士の積 V(st)*V(it) をプロットしていますが、両者のプロットが一致していることで乗算器が正しく動作していることが確認できました。

下のグラフは上のグラフと同じシミュレーション結果を対数軸で表示したものです。
対数軸でも V(stxst) と V(st)*V(it) のプロットは一致しており、乗算器の動作の正しさが追認されました。
V(it) がピークアウトした後、ほぼ直線的(指数関数的)に減少しています。

ここで V(st) が V(st) \( \ll 1 \) に収束するとき \( \gamma - \beta S(t) =k \) (=const) , \( I(t) = y \) とおくと

\( \frac{ \mathrm{ d } y }{ \mathrm{ d } t } =-ky \)

変数分離形の微分方程式として

\( \int \frac{ \mathrm{ d } y }{y} = -k \int \mathrm{ d } t \)

\( \log{|y|} = -kt + C_1 \)

これは \( I(t) \) が片対数グラフにおいて切片 \( C_1 \) (積分定数)、傾き \( -k \) の右肩下がりな直線で表される指数関数

\( y=C_2 e^{-kt} \)

に漸近することを示しています。


東京都コロナ第5波PCR陽性者数推移

下のグラフでは東京都のコロナ第5波のPCR陽性者数推移を曜日変動をキャンセルするためにPCR陽性者を7日間の移動平均で対数軸上に表示しています。

9月初めから10月末にかけて、50日間で1/100程度の割合で指数関数的に減少していることが判ります。

ここで \( t_1 \lt t_2 \) とすると

\( \log{ \left( \frac{I(t_1)}{I(t_2)} \right) } =-k(t_1 -t_2) \)

\( \log{100} = 50k \)

\( k= \frac{2 \log{10}}{50} = 0.092103 \)

参考文献*1の \( \gamma =0.1 \) という設定値の妥当性が、シミュレーションとこの計算結果からも確認できたのではないでしょうか。

以下、話は遡りますが、シミュレーション回路の積分器のコンデンサC2,C3,C4の容量の設定について説明します。
V(st) が 0 に収束した場合を考えると \( k= \gamma = 0.1 \) また \( \frac{ I(t_1) } {I(t_2 )} = 100 \) とすると

\( k= \frac{2 \log{10}}{t_2-t_1} =0.1 \)

\( t_2-t_1= \frac{2 \log{10}}{k} = 20 \log{10} =46.052 \)

タイムスケール1day→1msとするために、下のシミュレーション回路 ( \( \gamma =0.1 \) ) でのフィードバックによる \( I(t) \) のステップ応答が、46msで1/100になるようC3の値のカット＆トライを試みました。

下のシミュレーションのようにカット＆トライの結果、積分器のコンデンサC2,C3,C4の容量を995µFに決めました。


PCR新規陽性者の下げ止まり

東京都のPCR新規陽性者数は10月末までは指数関数的に減少していましたが、11月に入り20数人程度の日が続き下げ止まりとの声が聞かれるようになりました。
下のグラフはLTspiceによるSIRモデルシミュレーション回路のV3の感染爆発トリガのpulseオプションを外し、定常的に100µを入力するように変更したときのシミュレーション結果です。下げ止まりのような結果が現れています。

ここで下げ止まり後も \( \beta S(t) - \gamma \) が大きく変化していないという仮定のもと、東京都の外(国内外)から一日あたり何人程度の感染者が流入しているかを等比数列という観点から推定してみます。

項比 \( r ( 0 \lt r \lt 1) \) の等比数列第\(n\)項までの和を \( S_n \) とすると

\( S_n = a_0 \frac{1-r^n}{1-r} \)

\( \displaystyle \lim_{ n \to \infty } S_n = \frac{a_0}{1-r} \)

項比は、50日間で1/100に減少することから

\( r^{50} = \frac{1}{100} \)

\( r= \left( \frac{1}{100} \right)^{ \frac{1}{50}} =0.91201 \)

\( \displaystyle \lim_{ n \to \infty } S_n = \frac{a_0}{1-r} = \frac{a_0}{1-0.91201} = 11.365 a_0 \)

もしかすると東京都内へ今も感染者が平均2人/日程度流入し続けているのかもしれません。