数値シミュレーションによる検討

2020年3月21日版

文責:佐藤彰洋(ahsato@yokohama-cu.ac.jp)

遅れ付SIRモデルによるCOVID-19のパラメータ推定を各種データやパラーメータを変更して実行することで、新型肺炎の感染拡大抑止に関する方法について研究、分析、検討を行っていく。

数値シミュレーション結果のまとめ

  • 2020年2月28日頃から開始した大規模イベントの自粛の効果はそれ以前の感染率をベースラインとした場合、感染率を50%まで低下させていると見積もられる。
  • それに続く、2020年3月3日-4日から開始した小中高等学校の一斉休校の影響が生じていると考えられる2020年3月17日以降の感染者数を含んだ推計値はベースラインの59%まで感染率低減の効果と推計される。
  • 3月3日-4前後に海外で感染した感染者が帰国後に3月16日頃から多数国内の様々な地域で発症確認されているため、その増加分が小中高等学校の一斉休校措置で見込まれた減少分を打ち消していることがその原因ではないかと推察する。想定していなかった外部的要因が強く影響している可能性があるため、この施策の影響確認にはより長期間のデータを必要とし、引き続き確認のための追加推計を行わなければならない。
  • 現在実施されている感染率低減措置下における、感染率で、検査頻度を現在の6倍以上に高め、現在の6倍以上の頻度で環境中に存在する感染者を見つけ出し、隔離措置をおこなうことにより、感染拡大を終息させることが可能と予想される。
  • これまでの経済社会的活動で確認される感染率の0.07倍(7%)以下の感染率の値になるような大規模な感染率低減効果のある措置を広い領域で行うことで、感染者数を低減させ、感染拡大を終息させることができる(詳しくは社会的距離戦略における減少目標値の解釈を参照のこと)
  • 感染率を0.07倍(7%)以下とするような対策を14日以上継続して行うことで、感染者数の減少を確認することができるようになる
  • 感染者数の減少を確認してから感染者がいなくなるような終息状況となるのに必要とされる期間は、感染率を0.07倍(7%)以下に低減できるような大規模な対策を早めに行うほど短期間化される

遅れ付SIRモデルと変数の定義

(1)   \begin{eqnarray*}\frac{dS}{dt} &=& -\alpha(t)  \frac{S(t)I(t)}{N(t)} \\ \frac{dI}{dt} &=& \alpha(t-\tau)  \frac{S(t-\tau)I(t-\tau)}{N(t-\tau)} - \beta(t) I(t) \\ \frac{dR}{dt}  &=& \beta(t) I(t) \\ \frac{dD}{dt} &=& fatalityrate \times \beta(t) I(t) \\ \frac{dR_r}{dt} &=& (1-fatalityrate) \times \beta(t) I(t)  \\ N(t)  &=& S(t)+I(t)+R(t)  \end{eqnarray*}

(2)   \begin{eqnarray*} S(0) &=& N_0 \\ I(t) &=&  I_0(t) \quad (-\tau \leq t \leq 0) \\ R(0) &=& D_0 + R_{r0} \\ D(0) &=& D_0 \\ R_r(0) &=& R_{r0} \end{eqnarray*}

ここで、

  • S(t): 未感染の状態の人数
  • I(t): 感染した状態の人数
  • R(t): 回復または死亡した状態 の人数
  • D(t): 死亡者数
  • Rr(t): 回復者数

である。 感染率\alpha(t)と回復率(または除去率)\beta(t)は対策によりその値を変化させることができる。例えば、イベントを取りやめることや、直接接触の頻度を低下させて感染率\alpha(t)を段階的に減少させることができる。また、感染者を早急に見つけることで未感染者との接触を断つことにより、感染させる危険性を低下させ回復した人と同じ効果を得ることによりが可能である。この場合は\beta(t)を対策により上昇させることができる。fatalityrateは死亡率である。

今回のCOVID-19では、潜伏期間(3日~14日), 症状発症(1週間), 重症化(発症後2週間)を経てから感染が疑われて検査機会を獲得して、感染が確認できる。そのため、接触感染や飛沫感染などで、ウイルスを体内に取り込んで、実際に症状と確認できるまでの期間である潜伏期間(平均5.6日)と、発症後、風邪用の症状が確認される初期症状期間(7日)で、並びに、その後、特に重症化している場合に肺炎用の症状を示して、はじめてPCR検査を受ける機会が得られて、検査により陽性反応が出た場合に、感染者として確認できる状況となっている。

そのため、接触により罹患してから感染者として確認されるまでの遅れ時間τは14日~21日とするのが妥当(接触事象が発生してから実際に罹患して感染が確認できるようになるまで2週間以上必要)である。 遅れ時間τについては、12.6日~30日の間の値で、実際には数値シミュレーションを何通りかで実行することが望ましい。

パラメータの推定方法

各種の外部的に与えるべきパラメータは以下のように仮定する。

感染者数の日次時系列データI_o(t)が得られたとき、モデルから計算される数値シミュレーションにおける感染者数の値\hat{I}(t; \alpha, \beta)との間で二乗誤差を以下で定義する。 ここで、t_st_eはパラメータの適合を確認するための時間期間である。

(3)   \begin{equation*}E(\alpha, \beta) = \sum_{t=t_s}^{t_e} \Bigl(I_o(t)-\hat{I}(t; \alpha, \beta)\Bigr)^2\end{equation*}

パラメータ推定結果

感染率αと回復率βとを実際のデータからフィッテングにより推計してみる。感染者数の推移のデータは国内感染状況(NHK)まとめから取得した。

https://www3.nhk.or.jp/news/special/coronavirus/

日本における2月16日から3月14日までの感染者数のデータを使い、3月1日以前の感染者数をI_0(t)として、感染者数の初期値として使い、3月1日から3月14日までの感染者数の日次時系列データI_o(t)とモデル計算値\hat{I}(t; \alpha, \beta)との二乗誤差E(\alpha, \beta)をできるだけ小さくできるように 感染率\alphaと回復率\betaの推定を行った。 また、3月1日の回復者数と、死亡者数から、R_r(0)=31, D(0)=5, R(0)=36 を仮定した。

ここで、日次感染者数は積算での感染者数を用いた。これは、今回の新型コロナウイルスの感染者が回復と判断される条件(2回以上のPCR検査陰性)を満足しても、その後症状の再発がみられるため、暫定的条件での回復判断の数値を引くよりは、一度感染した感染者は今回のシミュレーション期間である数カ月は感染能力を有するというシナリオでシミュレーションを行うことで、もっとも最悪事象を想定した場合のシミュレーション結果を得るためである。最悪事象において確認できる各種の政策効果を判定することにより、それよりも緩和された回復者がある一定の割合で再発するという状況においては、それよりも政策効果が得られるため、より安全側の立場に立った判断をおこなうことができることは当然である。

二乗誤差ができる限り小さくなる\alpha\betaの値から、パラメータ推定値として以下の値を得た。

\alpha\betaE(\alpha,\beta)
0.4880.045523008.22
〇はモデルの数値シミュレーションから計算される感染者数、●は日次感染者数の実測値

致死率に関する検討

2020年3月1日現在の死亡者数は5名、2020年3月14日現在の死亡者数は21名となっている。この数値に一致するようにD(t)=fatalityrate \times R(t)とし、3月1日のD(t)の値が5名、3月14日のD(t)の値が21名と一致するようにfatalityrateを決定したところ、

fatalityrate=0.06

との値を得た。すなわち、モデル推計から導き出される死亡率は6%となる。

致死率がWHO報告の3.4%より高い理由としては、高齢者ほど重症化しやすいため、検査対象となりやすいことで、全体的な感染者数として認定される年齢階級が高齢者に偏りをしめしているため、死亡率を押し上げる原因となっているのではないかと考える。

この死亡率の推計により、死亡者数と回復者数の時間的な予測値を計算することができる。以下に死亡者数と回復者数の時間的発展を示す。

数値シミュレーションから得られた死者数の時系列(α=0.5, β=0.046, fatalityrate=0.06)
数値シミュレーションから得られた回復者数の時系列(α=0.5, β=0.046, fatalityrate=0.06)

長期的な傾向

上記で推計したパラメータを用いて、より長期的な傾向について、数値シミュレーションを行った結果を以下に示す。感染率低減抑制を行わず、このままの活動を継続し続けた場合は、感染者数I(t)、死亡者数D(t)、回復者数Rr(t)ともに単調に増加傾向を続けることが分かる。その増加速度は、感染者数、回復者数、死亡者数ともに、約2週間でおおよそ10倍の増加率である。

感染者数I(t)の時間発展の様子(y軸は対数スケールで表示)
死亡者数D(t)の時間発展の様子(y軸は対数で表示)

議論

データ期間の違いによるパラメータ推計値の変化

3月8日までの感染者はその14日前までである2月23日以前の経済社会的な活動の中で罹患者との接触事象が発生し罹患した人々である。この当時の活動は大型イベントが通常どおり行われており、企業活動も通常であった。

その後、2月26日から日本全国においてイベントの自粛や大型施設の閉鎖が始まったため、3月11日前後から感染率は減少していると期待がされ、その減少の効果によって実際に確認される感染者数は減少していると思われる。すなわち、感染率が低下しているはずである。では、実際2月26日から開始した、イベントの自粛や大型施設の閉鎖の対策には、どの程度の感染率の低下がもたらされているのであろうか。

これを確認するためにパラメータ推定に用いるデータ期間を3月1日から3月8日までと3月1日から3月13日までの2パターンで推計を行ってみる。

推定データ期間感染率 α回復率(除去率) β
2020年3月1日~3月8日0.6230.0505
2020年3月1日~3月14日0.5000.046

感染率は3月9日~3月14日までのデータを含んで推定したほうが、2月25日頃からのイベント自粛の影響が作用しているため、感染率の推定値が小さい。感染率は0.5 \div 0.623 \times 100= 80.2%の値に減少している。

更に、2020年3月12日頃から2020年2月29日~3月1日頃またはそれ以前に感染者と接触したことにより発症したともられる店舗従業員の発症事例が愛知県と北海道で複数確認されている。これは、それ以前に確認がされていた医療従事者とは異なる形でのこれまで想定されていなかった感染形態である。この種の感染形態は大規模イベントの自粛要請や、大規模公共施設の閉鎖措置では感染を未然に防ぐことが不可能であった感染事例である。

もう少しデータ期間を延長して確認していく必要があるが、3月14日時点で、イベント自粛の効果は直接接触の頻度を減少させ、感染率を約20%ほど低減される効果があったとみられる。他方、回復率(または除去率)βは若干悪化してきており、109.8%の値への悪化がみられる。これは感染者数が増加したことにより徐々に隔離作業に資源が大量に投入が必要となったことによる影響ではないかと推察する。

区分的感染率変動モデルによる施策の影響評価

さらに期間を延長し、2020年3月1日~2020年3月19日までのNHKまとめ積算日次感染者数(JPN-COVID-19-20200319.csv;2020年3月19日21時発表後にNHK報告値に更新があり2020年3月20日10時現在の値に修正した)を用いた区分的感染率変動モデルによるシミュレーション結果のフィッティングを行った。

外部的に与えられるパラメータのうち特に基本生産数R0については、ドイツで推計される最新の値を標準的な人間活動で生じるウイルス感染の値として仮定する。

区分的感染率変動モデルとは感染率を以下のように3つの期間に分割してそれぞれの期間とパラメータ推計値がベース推計値\alpha_0からどの程度の比率とするとシミュレーション結果と実測値の2乗誤差が小さくなるかによって、区間と感染率の減少程度を見積もる方法である。t_1=2020年2月28日(大規模イベント自粛要請後)、t_2=2020年3月5日(小中高等学校の一斉休校後)と仮定して\alpha_0, q_1, q_2を二乗誤差が小さくなるように選んだ。

(4)   \begin{equation*}\alpha(t) = \left\{ \begin{array}{cc} \alpha_0 & (t < t_0) \\ q_1 \alpha_0 & (t_0 \leq t \leq t_1) \\ q_2 \alpha_0 & (t_1 \leq t) \end{array}\right.\end{equation*}

パラメータ期間推定値
\alpha_0~2020年2月27日0.493
q_12020年2月28日~2020年3月4日0.5
q_22020年3月5日~0.59
数値シミュレーションから得られた感染者数〇と日次感染者数の実測値●

イベント自粛による感染率の低減効果はベースライン感染率を当初の50%まで減少させる効果はあったとみられる。また、小中高等学校の一斉休校実施後はベースライン感染率から59%程度の感染率低減効果と見積もられる。このことから2月末から行ってきた大規模イベントの自粛や公共施設の閉鎖、小中高等学校の一斉休校措置は感染率低減に効果を示していると考えられる。小中高等学校の一斉休校実施による感染率の低減効果があまり確認できていない背景として、ちょうどこの期間に外国旅行に出かけ、海外で感染した帰国者が帰国後に3月16日以降発症するという事例が多発していることが原因ではないかと想像する。この海外からの帰国者の発症事例の多発により、2020年3月3日以降で実施した小中高等学校の一斉休校の効果で国内で減少できた感染者数増加分以上の感染者数の増加が流入によりもたらされたと考えられる。

引き続き感染率低減効果のある直接接触頻度を低下させる施策を継続することで感染者数の増加速度を遅らせ、オーバーシュートの発生を抑制し、医療機関への負担を低減させながらクラスタ対策に資源を集中することで除去率(または回復率)\betaを上昇させる措置をとれる余裕を作りだし、終息を目指す戦略が有効と判断する。

日次感染者数の更新に伴って、パラメータ推定結果の変化の傾向を確認していくことにより感染率の低減効果がどの程度であったかについて継続して見積もりするめていく必要がある。

検査と隔離の頻度を高めた場合のシナリオ分析

検査を現状よりも広範囲かつ大規模に行うことにより環境中に存在する感染者を見つけ出し、隔離することにより他者への感染を抑制する努力目標について検討する。このような積極的な隔離措置は、除去率(または回復率)\betaの値を現状より強制的に大きくすることに対応する。

上記でおこなった2月末からの対策を織り込んだフィッティングより得られた、3区分的感染率を使って2020年3月20日から、検査頻度と隔離措置を大規模に行った場合に現状の\betaの値を何倍とすると感染者数の減少が確認されるようになるかについて数値シミュレーションにより検討する。

現状の検査と隔離の頻度を4倍に高めた場合
現状の検査と隔離の頻度を5倍に高めた場合
現状の検査と隔離の頻度を6倍に高めた場合

数値シミュレーションから、2020年3月20日から検査頻度を現在の6倍以上に高め、環境中に存在する感染者を現状より6倍多く探し出して隔離を行い、感染者が環境中でウイルスを未罹患者へ感染させない措置をとることができたとすると、感染者数の減少を認めることができるようになると予想する。

検査頻度と感染者の隔離頻度が現状の6倍以上のより大きな値であればあるほどに、早く終息をさせることが可能である。

すなわち、現状の感染率低減効果が認められる大規模イベントの自粛、公共施設の閉館、小中高等学校の一斉休校と並行して、検査頻度を高め、環境中に存在する感染者をより多く発見し(少なくとも現状の6倍以上の頻度で発見する)、隔離措置につとめることにより、新型コロナウイルスの感染拡大を終息させることが可能と予想する。

感染拡大低減のシナリオ分析

今後、更なる対策により、 α(t) = α0 Vul(t) : 感染率(ここで、 α0 は前述の感染率の推定値である) として、感染率をVul(t)を1以下に変更することで大幅な低減策を実施するシナリオを考える。この時、現状の感染率からどの程度の感染率減少をどの時点で行うのがよいかが問題となる。

(5)   \begin{equation*} Vul(t) = \left\{ \begin{array}{cc} 1 & (t< T) \\ q & (t \geq T) \end{array}\right.\end{equation*}

すなわちTqの値をどのように選ぶのが適切であるかという問題である。

どこまで感染率を減少させればよいか

このとき、Vul(t)をどの程度まで小さくできれば、感染者数は減少に転じるかをパラメータqを導入して検討してみる。仮定として3月15日にqをどの程度まで低減する措置を行いそれをその後継続すればよいかについてシミュレーションを行った。

T=3月15日に、q=0.2 (無対策から20%まで感染率を減少できた場合)と設定できる対策を実施した場合のシナリオでは、 その14日後の3月29日頃から感染者数の一旦の減少は見られるものの、その後増加傾向へ転じ感染拡大を抑止する効果はない。

T=3月15日に、q=0.1 (無対策から10%まで感染率を減少できた場合)と設定できる対策を実施した場合のシナリオでは、 感染者数の一旦の減少は見られるものの感染拡大を抑止することはできない。

T=3月15日に、q=0.08(無対策から8%まで感染率を減少できた場合)と設定できる対策を実施した場合のシナリオでは、感染者数増加に対する抑制効果に低減効果がみられる。

T=3月15日に、q=0.07(無対策から7%まで感染率を減少できた場合)と設定できる対策を実施した場合のシナリオでは、感染者数増加傾向を減少傾向に転じることが可能である。

T=3月15日に、q=0.06(無対策から6%まで感染率を減少できた場合)と設定できる対策を実施した場合のシナリオでは、感染者数増加傾向を顕著に減少傾向に転じることが可能である。

これらの分析の結果から、都市の封鎖、多くの人々の移動制限、並びに、自宅待機をほぼ同時に行うようなある程度強制力のある、大幅な効果が期待できる感染率低減措置をとらない限りは感染者数の減少傾向へ転じるような効果が期待できないと結論付られる。社会的距離戦略を採用する場合は、現在人々が行っている直接接触頻度を現状の6%-7%程度まで低減させることを2週間以上に渡り継続しなければならない。

また、前述の検査頻度を最大限に高め、感染者を発見し、隔離措置を行うことで、感染者と未罹患者との直接接触を限りなく減少させる措置も並行して行うことで、感染率を低下させるとともに、除去率(または回復率)\betaを増加させることで、感染者数の減少を認めることが容易になると予想する。

感染率低減対策の導入時期と終息までに要する期間との関係

次に、対策日Tの効果について次に検討する。q=0.07と固定して、対策日Tをそれぞれ3月15日、3月19日、3月25日と変化させた場合のシミュレーションの結果を以下に示す。

3月19日から効果的対策を講じた場合、3月15日に対策を講じるのと比較して、ピークでの感染者数が大きくなり、ピークアウトする時期が4日遅れることになる。

3月25日から効果的対策を講じた場合、3月15日に対策を講じるのと比較して、更にピークでの感染者数が大きくなり、ピークアウトする時期が10日遅れる。ピークアウトする感染者数が大きくなり、ピークアウトする時期も遅れるためその後終息するまでに措置を継続しなければならない期間が長期化ことは容易に理解できる。

上記の状況を定量的に理解するために、感染者数が減少に転じてから、感染者数1000名を下回るまでの期間を、対策実行日の関数として確認する。以下表では、q=0.07の場合に、対策実行日と終息確認に要する日数を示している。

対策実行日ピーク確認日ピーク感染者数(人)終息確認日日数(日)
2020-03-152020-03-2832292020-11-13243
2020-03-192020-04-0146972021-2-12339
2020-03-252020-04-0780952021-6-27446

そのため、できるだけ感染者数が少ない人数が確認された時点で、早めに規模の大きな対策を行い、感染率をこれまでの経済社会的な活動で確認される値の0.07倍(7%)以下まで減少させるような対策を講じる方が、終息までに要する経済社会的なコスト(人的資源、金銭的費用)は小さくできる。

次にq=0.05の場合の対策実施日と終息確認日までの期間を調べてみる。以下表に対策実施日ごとのピーク感染者数と対策を実施してから終息が確認されるまでの日数を示す。q=0.07よりもq=0.05のほうが対策実施から終息が確認されるまでの日数が短縮化されていることが分かる。また、ピーク感染者数もわずかながら減少することが確認できる。

対策実施日ピーク確認日ピーク感染者数(人)終息確認日日数(日)
2020-03-152020-03-2832232020-06-1289
2020-03-192020-04-0146892020-07-14117
2020-03-252020-04-0780812020-08-31159

このとこから、感染率を低減させる効果のある対策をできるだけ厳密にすることにより(qの値をより小さくする)、終息するまでの時間は短縮されると判断される。