はじめに

　新型コロナウイルス感染症(COVID-19)の感染拡大を防止する取り組みとして、都市封鎖、行動・営業自粛、RT-PCR検査やクラスター対策、接触者追跡アプリ、感染者行動監視アプリ、などの対策が様々な国で取り組まれています。日本はこれまで、PCR検査、クラスター対策班による追跡、緊急事態宣言による行動・営業自粛を主軸に取り組んで来ましたが、2020年5月6日時点で、中国、イタリア、イラン、アメリカ等で見られた爆発的感染拡大には至っていません。しかし、このまま行動・営業自粛を継続し続けることは、経済・教育・文化など様々な側面で弊害があり、出口戦略を模索する段階に来ています。ワクチンの開発が積極的に行われていますが、安全性・効果の確認、生産体制の確立などを考慮すると、早くて１年以上の期間が必要と目されています。
　医療崩壊を防ぎつつ、段階的な行動・営業の制限解除を行うには、より効率的な感染者の隔離対策が必要です。その重要なアイテムとして、接触者追跡アプリを導入することによる効果をシミュレーションにより検証してみたいと思います。

世界各国の状況

　まず、世界各国の現在のCOVID-19の医療面での状況を見てみます。下図は、感染者が１万人以上の国において、横軸を回復率（回復者数／要治療者数）、縦軸を死亡率（死亡者数／要治療者数）とした平均値のグラフです。計算方法はこちらの記事参照。

　日本は回復率1.5％、死亡率0.25％付近の位置にいます。回復率で言えば、フランス、ベルギー、シンガポールと同等で、死亡率でいえば、ドイツ、オーストリア、パキスタン、ベラルーシと同等です。
　この中で、中国やイスラエル、シンガポール、韓国などのコロナ追跡アプリが成果を挙げているとされていますが、感染源とされ初動の遅れた中国を除いて、確かにイスラエル、シンガポールや韓国の死亡率は低い水準です。
　既に世界各国で１０種類以上のコロナ追跡アプリが開発されており、
（参考：Top 10 popular smartphone apps to track Covid-19)
　日本でもAppleとGoogleで共同開発したAPIを使ったアプリを開発中のようです
（参考：コンタクト・トレーシング・アプリの開発に関して）。
　AppleとGoogleは、コロナ追跡アプリのストア登録には各国の当局の承認が必要としていますので、複数アプリが乱立することはないと思います。
（参考：各国のコロナ対策アプリはどのようにプライバシー対策を行なっているのか？）

計算の前提

　さて、このように導入が進められているコロナ追跡アプリですが、期待される効果として、

濃厚接触者を素早く検知して、検査・隔離・治療を早期に促す効果（特に無症状感染者）
聞き取り調査の際に「感染経路不明」となるケースの低減
クラスターの早期検知
感染者をリスクの高い場所に入場させないための通行証代わり
その他、疫学的に有益な知見の蓄積など

など様々な効果が期待されますが、今回は特に１に注目して効果を定量的に評価したいと思います。

　その際、計算モデルとして、以前の記事で示した、クラスター基準SEIRモデルを使います。このモデルを使う主な理由は次の通りです。

COVID-19の感染はクラスターが主体となって引き起こされているという知見に合うこと
実効再生産数$R_t$の時間変化を観測値から推定して織り込んでいること
SEIRモデルを現実に当てはめる際に課題となるSの初期値に依存しないこと

特に、通常のSEIRモデルでは基本再生産数$R_0＝2.5$の場合、ピークアウトが生じるのは感染者が人口の60%程度($1-1/2.5=0.6$)などと言われていますが（集団免疫の獲得）、１００万人当たりの感染者数が多い国を見ても、カタールやスペインで5-6千人（/100万人）ですので、**感染率はせいぜい0.6％**と見積もられます。**日本ではわずか0.01%です。**なので、観測されているCOVID-19のピークアウト現象を単純なSEIRモデルでは説明できないと考えるのが現時点では妥当と思われます。

計算モデル

　計算式を示します。まず、変数を導入します。

S:感染症に対して免疫を持たない者(Susceptible)
E:感染症が潜伏期間中の者(Exposed)
I:発症者(Infectious)
R1:アプリ以外のルートで感染症から回復し免疫を獲得した者(Recovered1)
R2:アプリ検知により隔離された者(Recovered2)
R:感染症から回復し免疫を獲得した者(Recovered)。($R=R_1+R_2$)

パラメータとして、次を導入します。

$\alpha$:アプリ隔離率。アプリによって隔離される１日あたりの潜伏期間中および発症中の感染者の比率と定義
Rp(t):時刻tにおける実効再生産数
lp:潜伏期間（5日を採用）
ip:感染期間（8日を採用）
Th:実効再生産数の半減期（7.5日を採用）

次式が、アプリによる隔離効果を組み込んだクラスター基準SEIRモデルです。

\begin{eqnarray}
\frac{dRp}{dt} &=& - \frac{ln2}{T_h}Rp \\
\frac{dS}{dt} &=& -\frac{1}{ip} Rp \cdot I \\
\frac{dE}{dt} &=& -\frac{dS}{dt} - \frac{1}{lp} E - \alpha E\\
\frac{dI}{dt} &=& \frac{1}{lp} E - \frac{1}{ip} I - \alpha I\\
\frac{dR_1}{dt} &=& \frac{1}{ip} I \\
\frac{dR_2}{dt} &=& \alpha E + \alpha I \\
\frac{dR}{dt} &=& \frac{dR_1}{dt} +\frac{dR_2}{dt} 
\end{eqnarray}

ちなみに、次式が成立しています。

\frac{d}{dt}(S+E+I+R_1+R_2)=\frac{d}{dt}(S+E+I+R)=0

ポイントは、アプリ効果によってEからもIからも一定比率で隔離されることです。

計算結果

　今回は、結果から先に示したいと思います。実効再生産数の初期値を$R_p(0)=10$とします。また、初期状態を$(E,I,R)=(0,1,0)$とします。

アプリ隔離率α=0のとき

アプリによる隔離が全くない場合です。最終的に113人が感染します。

アプリ隔離率α=0.1のとき

アプリによる隔離率が10[%/日]の場合です。つまり、１日あたり10％の感染者が隔離されると想定します。
最終的に26人が感染します（R1+R2）。感染者が77％減少しました。

アプリ隔離率α=0.2のとき

アプリによる隔離率が20[%/日]の場合です。
最終的に11人が感染します（R1+R2）。感染者が90％減少しました。

アプリ隔離率αと最終的な感染者数の関係

アプリによる隔離率$\alpha$と最終的な感染者数$R(T_{max}), R_1(T_{max}), R_2(T_{max})$の関係をグラフにしてみます。
およそ**$\alpha=0.07$以上で、アプリによる隔離者の方がアプリ以外の隔離者を上回って**いきます。

考察

以上から、クラスター基準のSEIRモデルに基づく、コロナ追跡アプリによる早期隔離の効果に関して、次の傾向がシミュレーションから導けます。

アプリによる隔離率が20[%/日]以上であれば、最終的な感染者数を90%以上削減する効果が期待できる。
アプリによる隔離率が7[%/日]以上であれば、アプリによる隔離者の方がアプリ以外の隔離者を上回る。

さらに言えば・・・

韓国の感染者数が低く抑えられている主要な要因は、PCR検査数よりもアプリ効果の方が大きい可能性があります。
同様に、日本の感染者が低く抑えられている主要な要因は、クラスター対策班や保健所によるクラスター追跡の効果である可能性が高いです。国民性などの他の要因も十分考えられますが、日本は世界でも屈指の高人口密度国であるという不利な条件も十分加味する必要があります。
PCR検査の際の人的（医師、看護師、検査技師、保健所職員etc.）、物的資源（防護服、マスク、検査キットetc.）が有限であることを考慮すると、抗体検査や抗原検査など非接触で簡易かつ高感度・高特異度な検査の開発が望まれます。コロナ追跡アプリは仮に検査しなくても、特に無症状者を隔離できる点で有益です。
コロナ追跡アプリは、本記事で示した通り絶大な効果が期待できますので、プライバシーに十分配慮しつつ積極的に導入すべきだと思います。イベントや施設への通行証機能も持たせればクラスターの予防も期待でき、経済活動再開へのキーアイテムともなりえると思います。

参考リンク

下記のページを参考にさせて頂きました。

Top 10 popular smartphone apps to track Covid-19
コンタクト・トレーシング・アプリの開発に関して
各国のコロナ対策アプリはどのようにプライバシー対策を行なっているのか？
新型コロナウイルスの医療崩壊仮説を検証してみる
新型コロナウイルスの日本と他国との広がり方の違いをクラスター対策モデルで説明してみる

Pythonコード

最後に、本記事で使用したPythonコードを貼っておきます。

import numpy as np
import matplotlib.pyplot as plt

# ODE solver
def my_odeint(deq, ini_state, tseq, keys):
    sim = None
    v = np.array(ini_state).astype(np.float64)
    dt = (tseq[1] - tseq[0])*1.0
    for t in tseq:
        dv = deq(v,t, keys)
        v = v + np.array(dv) * dt
        if sim is None:
            sim = v
        else:
            sim = np.vstack((sim, v))
    return sim

#define differencial equation of seir model
def seir_eq8(v, t, keys):
    Th = keys['Th']
    lp = keys['lp']
    ip = keys['ip']
    al = keys['al']
    #
    Rp  = v[0];
    s   = v[1];
    e   = v[2];
    i   = v[3];
    r1  = v[4];
    r2  = v[5];
    r   = v[6];
    #
    dRp  = - np.log(2)/Th * Rp
    ds   = - Rp/ip * i
    de   = - ds     - (1/lp) * e - al * e
    di   = (1/lp)*e - (1/ip) * i - al * i
    dr1  = (1/ip)*i
    dr2  = al * e + al * i
    dr   = dr1 + dr2
    return [dRp, ds, de, di, dr1, dr2, dr]

# simulation 1
def calcsim(Rp0, keys):
    # solve seir model
    #          Rp  S  E  I  R1 R2  R
    ini_state=[Rp0, 0, 0, 1, 0, 0, 0]
    t_max=180
    dt=0.01
    t=np.arange(0,t_max,dt)
    #
    sim = my_odeint(seir_eq8, ini_state, t, keys)
    #
    plt.rcParams["font.size"] = 12
    fig, ax = plt.subplots(figsize=(10,5))
    ax.plot(t,sim[:,[2,3,4,5]]) # extract E I R1 R2
    ax.set_xticks(np.linspace(0,t_max,19))
    yw = 10; yn = int(120/yw)+1
    ax.set_yticks(np.linspace(0,yw*(yn-1), yn))
    ax.grid(which='both')
    ax.legend([ 'Exposed', 'Infected','Recoverd1(normal)', 'Recoverd2(App-based)'])
    ax.set_xlabel('date')
    ax.set_ylim(0,)
    plt.show()
    print("R(tmax):{}".format(sim[-1,6]))

# do simulation 1
keys = {'lp':5, 'ip':8, 'Th':7.5, 'al':0.0 }
calcsim(10, keys)
keys = {'lp':5, 'ip':8, 'Th':7.5, 'al':0.1 }
calcsim(10, keys)
keys = {'lp':5, 'ip':8, 'Th':7.5, 'al':0.2 }
calcsim(10, keys)


# simulation 2
def calcAltoRfin(Rp0):
    #          Rp  S  E  I  R1 R2  R
    ini_state=[Rp0, 0, 0, 1, 0, 0, 0]
    t_max=180
    dt=0.01
    t=np.arange(0,t_max,dt)
    #
    lp = 5
    ip = 8
    keys = {'lp':lp, 'ip':ip, 'Th':7.5, 'al':(0.) }
    #
    mklist = lambda e, l : np.append(np.array(e),l)
    #
    rslt = []
    for i in np.linspace(0,0.4,20):
        keys['al'] = i
        sim = my_odeint(seir_eq8, ini_state, t, keys)
        r = mklist(keys['al'], sim[-1][4:])  # (i, R(tmax), R1(tmax), R2(tmax))
        rslt.append(r)
    #
    rslt = np.array(rslt)
    ymax = max(rslt[:,2])
    plt.rcParams["font.size"] = 12
    fig, ax = plt.subplots(figsize=(10,5))
    #
    ax.plot( rslt[:,0], rslt[:,3], 'b')
    ax.plot( rslt[:,0], rslt[:,1], 'g')
    ax.plot( rslt[:,0], rslt[:,2], 'r')
    ax.legend([ 'R(Total infected)', 'R1(normal purge)', 'R2(App-based purge)'], loc='upper right')
    ax.grid(which='both')
    ax.set_xlabel('application-based purge rate [/day]')
    ax.set_ylabel('total infected cases at tmax')
    yw = 10; yn = int(120/yw)+1
    ax.set_yticks(np.linspace(0,yw*(yn-1), yn))
    #
    for tat in [0,0.1,0.2,0.3,0.4]:
        idx = [i for i in range(len(rslt[:,0])) if rslt[i,0] >= tat][0]
        print("R_fin with Al:{} is {}".format(rslt[idx,0], rslt[idx,3]))
    #
    ax.set_xlim(0,)
    ax.set_ylim(0,)
    plt.show()

# do simulation 2
calcAltoRfin(10)