病原体ゲノムデータから抗菌剤耐性のフィットネスコストとベネフィットを推定する
病原体ゲノムデータから抗菌剤耐性のフィットネスコストとベネフィットを推定する
https://royalsocietypublishing.org/doi/10.1098/rsif.2023.0074
デイビッド・ヘレカル
,
マット・キーリング
,
ヨナタン・H・グラッド
および
グザヴィエ・ディドロ
発行：2023年6月14日https://doi.org/10.1098/rsif.2023.0074
要旨
多くの細菌性病原体集団における抗生物質耐性の増加は、公衆衛生に対する大きな脅威である。抗生物質に対する耐性は、細菌がこの抗生物質にさらされたときにフィットネス上の利益をもたらすが、耐性はまた、感受性の相手と比較して耐性病原体にコストをもたらすことが多い。私たちは、多くの細菌病原体や抗生物質について、耐性による利益とコストを十分に理解していませんが、それらを推定することで、耐性の拡大を抑制または防止する方法で抗生物質をよりよく使用することにつながります。ここでは、耐性菌のコストとベネフィットを明示的にパラメータとして含む、感受性変異体と耐性変異体の共同疫学に関する新しいモデルを提案する。このモデルでは、感受性系統と耐性系統の系統データを用いてベイズ推論を行うことができ、両者のデータを組み合わせることで、耐性コストと耐性ベネフィットのパラメータを分離して個別に推定できることを示す。我々は、推論手法をいくつかのシミュレーションデータセットに適用し、優れたスケーラビリティと精度を実証した。米国で2000年から2013年にかけて収集された淋菌ゲノムのデータセットを分析した。その結果、フルオロキノロンに耐性を持つ無関係な2つの系統が、類似した流行ダイナミクスと耐性パラメータを共有していることがわかった。フルオロキノロンは耐性レベルの上昇により淋病の治療が放棄されたが、今回の結果から、耐性が再び拡大することなく、10%程度の少数例の治療に使用できる可能性が示唆された。
はじめに
多くの病原体の抗菌薬耐性レベルは、過去数十年の間に心配なほど上昇している。CDC（米国疾病対策センター）が発表した抗菌薬耐性による脅威に関する報告書では、淋菌を含む3つの微生物が緊急の脅威レベルに分類され、さらに12の微生物が公衆衛生に対する深刻な脅威とされています [1]。抗菌薬耐性に関するレビューでは、耐性によって世界中で少なくとも年間70万人の命が奪われ、現在の傾向が続くと2050年までに死者数が年間1000万人に上る可能性があると推定されており［2］、最近の研究では、2019年には耐性に関連する死者が500万人近くいると推定されています［3］。1970年代以降、新しい抗菌薬が開発・配備されたことはほとんどなく、一方、新薬に対する耐性は最初の導入後すぐに出現することが多いため[4]、いくつかの病原体は完全に治療不能になる危険な状態に陥っています。抗菌薬耐性に効果的に対処するには、耐性出現の原因となる疫学的・進化的要因や、病原体集団における耐性の拡散をより深く理解する必要があります。この目標を達成するためには、抗菌薬耐性に関する数理モデルを開発し、情報量の多い観察結果を用いて疫学モデルをしっかりと統計解析する必要があります。耐性に関するこのようなモデリングアプローチは1990年代後半に開始され[5,6]、異なる生物、拡散様式、研究規模や状況に適した多くのモデルの開発につながった[7]。
耐性は、抗菌薬の存在下でそれを獲得する病原体に明確な体力的利益をもたらす。したがって、この体力的利益の純価値は、特定の抗菌薬が使用される頻度が高いほど、病原体そのものに対して、あるいは無症状で運ばれる病原体の場合、より一般的に増加する。しかし、耐性は通常、病原体にとってコストとなる。この効果を最も簡単に示すのは、抗菌薬の使用を中止することで耐性率が低下する場合である。しかし、多くの病原体や抗菌薬において、耐性がもたらすコストと利益はまだ十分に理解されていない[9]。耐性菌の蔓延を食い止め、あるいは逆転させることを期待して適性コストを利用するために提案された公衆衛生介入策の潜在的有効性を評価するための確かな根拠を提供するためには、耐性菌の利益とコストのより良い定量化が必要である [9]．例えば、最近、イングランドにおける10年間のセフィキシムに対する感受性および耐性の淋病患者数が、この抗生物質に対する耐性に関連するコストと利益を定量化するために分析された［10］。これらの推定値を用いて、セフィキシム耐性レベルの上昇を招くことなく、淋病患者の少数派（約25％25％）の治療にセフィキシムを再導入でき、現在使用されている抗生物質に対する耐性出現のリスクを低減できると予測された。さらに、耐性の適性コストの程度はゲノム背景によって異なる可能性があり[11]、耐性の適性コストを利用しようとする介入の効果は系統に依存する可能性があるようなものである。そのため、系統ごとのフィットネスコストを推定する必要がある。本研究の目的は、処方方針の変更が特定の耐性系統の集団動態に与える寄与を定量化することである。これは、耐性の全体的な生態や耐性表現型の最終的な運命に興味を持つ研究（例えば[12]）とは対照的である。
病原体ゲノムデータは、感染症の進化的・疫学的ダイナミクスの理解に役立つ大きな可能性を秘めている[13]。この系統力学的アプローチの重要な利点は、ゲノムデータの解析がサンプリングバイアスの影響を受けにくいことで、特にサンプリングを条件として祖先のプロセスを記述するcoalescentフレームワークを使用する場合[14]にある。いくつかの研究では、抗菌剤耐性に関連するフィットネスコストに光を当てるために、このアプローチを使用している。例えば、ある研究では、メチシリン耐性黄色ブドウ球菌の系統の成長率とβラクタム薬の消費量との関連性が示された [15] 。また、HIV [16]や結核菌 [17]における耐性変異の相対的な伝播適性を定量化した研究もある。ここでは、感受性系統と耐性系統の系統力学的軌跡を、一定である体力コストと抗菌薬消費量に依存する体力ベネフィットの関数として明示的にモデル化することによって、異なるアプローチをとることにする。つまり、抗菌薬の使用量、感受性系統のゲノムデータ、耐性系統のゲノムデータの3つを入力とする。そこから、耐性のコストと利益を分離し、系統力学的な軌道を予測するのに必要なパラメータを提供し、耐性の脅威を悪化させずに抗菌薬を使用する方法に関する推奨事項を提供します。全体として、私たちが興味を持っているシナリオは、大規模な集団レベルでの全体的な耐性ダイナミクスである。このようなシナリオでは、耐性菌の発生は輸入よりもむしろ地元での感染に起因することが多いでしょう。この論文で紹介した方法は、輸入品や複雑で異質な感染経路が支配的な小さな集団、例えば病院内の院内感染などに適用するつもりはありません。このようなシナリオでは、出生-死亡型モデルを用いた別のアプローチがより適切であろう [16,17] 。方法
2.1. 全体的なアプローチ
病原体の系統データには、研究対象の病原体の過去の個体数動態に関する情報が含まれている[13,18]。流行過程が単純なコンパートメント型流行モデルで十分に特徴付けられるという前提のもと、この集団サイズダイナミクスに関する情報を流行軌跡に変換することができる [19,20]．これらの疫学的軌道は、特定の抗菌薬に対する耐性のフィットネスコストとベネフィットの効果を考慮した疫学モデルを用いて記述することができる。この抗菌薬の使用が時間と共に変化すると、考慮する特定の系統の純適性も変化する。その結果、流行の軌道が変化することになる。しかし、流行軌跡のすべての変化が、耐性表現型のフィットネスの変化に起因するわけではありません。感受性の枯渇や宿主の行動の変化などの交絡要因も、流行の軌跡に影響を与えるだろう。以下に詳述する比較的穏やかな仮定の下では、これらの交絡因子の変化は他の系統にも等しく影響を与えることになる。そこで、理想的には近縁で、一次治療として大量に使用される他の抗菌薬に対する耐性プロファイルが同じである感受性系統のデータを「対照」として使用することができます。感受性の高い系統と耐性の高い系統の軌道の違いは、耐性に起因するものであり、関連するフィットネスコストとベネフィットのパラメータを推定することが可能である。
ある病原体が、ある抗菌剤で治療されている、あるいは治療されていた大規模な集団のレベルで感染症を引き起こしていると考えてみよう。過去のある時点で、この抗菌剤に耐性を持つ系統が1つまたは複数生まれたと仮定する。我々の目的は、ある系統がこの抗菌薬に耐性を持つことによる適合コストと利益を、時間経過に伴う対象抗菌薬の使用量の関数として定量化することである。そのためには、この病原体による感染症を治療するための、与えられた抗菌薬の経時的な使用を定量化するデータと、経時的にこの病原体による感染症から分離された配列決定済みの症例株の妥当なサンプルが必要です。さらに、個々の分離株の耐性プロファイルを特徴づける情報が必要であり、これはin vitroでの耐性スクリーニングによって得られるか、あるいはin silicoで配列から予測することができる[21]。これらのサンプルの系統樹は、例えばBEAST [22]、BEAST2 [23]、BactDating [24]などを用いて推定されます。この系統樹を解析の出発点として[25]、どのサンプルが耐性および感受性系統に属するかを特定し、目的の抗菌薬に完全に耐性または感受性であるが、その他の耐性プロファイルが類似している関連系統をさらなる研究のために選択する。なお、耐性モデリング研究[7]では、簡略化のため、耐性は2値形質として扱い、サンプルは抗菌薬に対して耐性または感受性のどちらかであることが通常である。
2.2. 感染モデルの導出
抗菌薬耐性のフィットネスコストとベネフィットを推定するためには、トランスミッションモデルを指定する必要がある。我々は、以前の感染によって再感染に対する免疫が与えられない場合に、ある治療抵抗性の表現型を持つ特定の系統の適合性パラメータを推定することに焦点を当てる。宿主集団が非構造的であり、過去の感染によって免疫が付与されないという単純化した仮定の下では、多系統の感受性-感染-感受性（SIS）が妥当なモデルである [26,27]. このモデルは、より一般的には多系統SISと呼ばれる。異なる系統の保有量の変動は、宿主の人口動態や行動の変化などの外的要因に起因することもあり得ます。このような変動が考慮されないまま放置されると、ある抗菌薬に対する耐性のコストとベネフィットの推定に偏りが生じる。そこで、感染率β(t)と集団サイズN(t)を時間的に変化させるモデルに変更することにした。これにより、以下のn連立常微分方程式(ODE)の系で記述されるn連鎖モデルへと発展する：
anddI1(t)dt=β(t)S(t)I1(t)N(t)−γ1(t)I1(t),dI2(t)dt=β(t)S(t)I2(t)N(t)−γ2(t)I2(t),⋮dIn(t)dt=β(t)S(t)In(t)N(t)−γn(t)In(t),⎫⎭⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪d�1(�)d�=�(�)�(�)�1(�)�(�)−�1(�)�1(�), d�2（�）d�＝�（�）�2（�）�（�）-�2（�），⋮d�（�）d�＝�（�）�（�）�（�）-�（�）�（�）， }．
2.1
ここで、Ij(t)は時刻tにおけるj番目の系統の感染者数を表し、β(t)はどの系統にも特有ではない変化、例えば宿主の行動などにより時間と共に変化する感染率である。N(t)は宿主の個体数で、これも人口動態的な要因で時間と共に変化する可能性がある。γj(t)は、時間tにおけるj番目の系統の回復率である。これらは、時間と共に変化する抗菌剤の使用量に依存するため、時間と共に変化する場合もしない場合もある。最後に、S(t)は感受性宿主の数を表す：
S(t)=(N(t)-∑ j=1nI j(t)).�(�)=(�(�)-∑ �=1�� �(�)).
2.2
一般的に、このモデルは、抵抗性プロファイルが表現型的に類似しているすべての系統を平均化する、単純に2系統モデルに縮小することができる。しかし、これは望ましくない。同じ抵抗性表現型を持つ系統の中には、ゲノム背景が異なるためにフィットネスが異なるものがあり、それが推定を混乱させるからである。さらに、この種のモデルはゲノムの枠組みでは扱いにくい。なぜなら、系統データは一般に特定の系統の動態についてのみ情報を提供することになるからである。このことは、解析が研究対象の系統に対して有効であり、ある病原体の抵抗性の全体的な動態に外挿することができないことも意味します。
したがって、我々は個々の系統の解像度に注目する必要がある。宿主の集団サイズや行動の変動などの環境効果は、集団がうまく混合されていれば、すべての系統に等しく影響することに留意する。これらの効果の組み合わせをb(t) = β(t)S(t)/N(t) とする。b(t)の知識軌跡を条件とすると、式(2.1)のODEは非結合となり、これによって、これから注目する系統に対応する非結合方程式に体系化することができる。そのため、b(t)は推論が必要なランダムなオブジェクトとして扱うことにする。さらに、感受性系統の平均回復率γsは時間と共に変化しないと仮定する。一方、耐性系統のそれは、ある患者が目的の抗菌薬で治療された場合はγT = qT + γs、それ以外の場合はγU = qU + γsの二つの値のどちらかを取る。また、登録された症例が目的の抗菌薬で治療された割合u(t)を考慮すると、耐性系統の平均回復率は次のように完全に決定される。
γr(t)=u(t)γT+(1-u(t))γU.��(�)=�(�)�+(1-�(�))��.
2.3
これで、感受性の高い系統と抵抗性の高い系統のそれぞれについて、これから使うモデルの方程式を完全に書き下すことができます：
dIs(t)dt= anddIr(t)dt= b(t)Is(t)−γsIs(t)b(t)Ir(t)−[u(t)γT+(1−u(t))γU]Ir(t).⎫⎭⎬⎪⎪⎪⎪⎪⎪d��(�)d�= �(�)��(�)−����(�)andd��(�)d�= �(�)��(�)−[�(�)��+(1−�(�))��]��(�).}
2.4
実際には、すべての症例に目的の抗菌薬を投与した場合と、目的の抗菌薬を全く使用しなかった場合の感受性系統と耐性系統の回復率の差に興味があります。これらを次のように表す。
とqT=γT-γsqU=γU-γs.}��=��-�� と��=��-��.} である。
2.5
したがって、qTはqT＜0の場合の抵抗の体力的利益を、qUはqU＞0の場合の抵抗の体力的コストをとらえていると解釈できる。
このモデルは、尤度に系統に関連した項を追加し、必要なパラメータを追加するだけで、任意の数の耐性系統と感受性系統に適用することができる。これは、個々の系統がb(t)の条件付きで独立しているので簡単であるが、簡単のために、残りの方法の説明は、単一の感受性系統と単一の耐性系統の場合に焦点を当て、一般的なケースは、簡単に拡張することができる。
2.3. 系統樹へのリンク
疫学モデルを定義したので、次はそれを系統発生プロセスにリンクさせることができる。19,28]に基づき、一対の系統の瞬時合体率を次のように導出することができる。
λs(t)=2b(t)Is(t)，λr(t)=2b(t)Ir(t)��(�)=2�(�)�(�)，�(�)=2�(�)�(�) �� (�)
2.6
を、それぞれ感受性集団と耐性集団に適用した。したがって、時刻s1＜-＞＜snでn枚の葉、時刻c1＜-＞＜cn-1でn-1枚の合体事象、時刻tでA（t）系統を持つ年代別系統樹gの尤度はGriffiths & Tavare［29］で与えられる：
p(g|λ(t))=exp(-∫∞-∞1A(t)≥2λ(t) dt)∏i=1n-1λ(ci), �(�|（�）)=exp(-∫∞-∞1[�（�）≥2](A（t）2）�（�） d�)∏=1�-1�（��)、
2.7
ここで、λ(t)=λs(t)、λ(t)=λr(t)はそれぞれ感受性系統、耐性系統の場合である。しかし、ほとんどの場合、そして実際に我々の場合、式(2.7)の積分は解析的に難解なものではありません。さらに、抗生物質の使用データが系統樹全体に及ぶことはまずない。そこで、抗生物質の使用データと研究対象の系統樹の交点である[tmin, tmax]で切り捨てた系統樹の近似尤度を定義します。
tmin、tmax]の区間を、ti - ti-1 < Δtとなるような細かいメッシュtmin = t1 < t2 < t3 < - - < tN = tmaxに分割し、tminからtmaxまでのすべてのサンプリングと合体時間をメッシュに含めるのである：
p(g|λ(t))=exp(-∑i=2N(ti-ti-1)(A(ti-1)2)λ(ti-1))∏i=1n-11[ci∈[tmin,tmax]] λ(ci). �(�) =exp(-∑�=2�(�-�-1) (�（�-1）2) �（�-1))∏�=1-11[� ∈［�min，�max］］�（��).
2.8
系統と疫病の関係をどのように扱うかというアプローチは、決定論的疫病モデルによって決定される移動と時間変化のないNe(t)を持つ構造化合体というのが事実であることに注意されたい。私たちと同じようなアプローチは、選択の有無にかかわらず、突然変異の予想年齢を正式に研究するために使用されています[31]。しかし、その場合、集団は異なる対立遺伝子に対応し、Ne(t)曲線は、Wright-Fisher拡散によって時間的に前方に決定される、与えられた対立遺伝子を持つ集団の割合に従う。また、個々の対立遺伝子に対応するデメ間の移動は、さらに組換えに対応して加えることができる[32]。
2.4. ベイズ推論
まず、時間を区間[tmin, tmax]から[ -1, 1]に再スケールする。この再スケーリングに関連するスケールファクターD = (tmax - tmin)/2 を、γs¯=γsD�¯=��� と定義してモデルで考慮することにする。
モデルは、時変係数を持つ各系統の独立した一次線形同質ODEで構成されている。初期条件Is(0)=Is0、Ir(0)=Ir0を前提とした時刻tでの解は、時刻tまでの瞬時レートの積分で求めることができる：
となり、Is(t)=Is0exp{∫t0b(τ)-γs dτ} となる。 Ir(t)=Ir0exp{∫t0b(τ)−[u(τ)γT+(1−u(τ))γU] dτ}.⎫⎭⎬⎪⎪⎪⎪⎪⎪⎪⎪��(�)=��0exp{∫0��(�)−�� d�} and��(�)=��0exp{∫0��(�)−[�(�)��+(1−�(�))��] d�}.}
2.9
このモデルは、b(t)の事前分布を適切に選択することが難しく、また初期条件とb(t)の間の依存構造が非常に複雑であるため、推論を行うには適していません。そのため、Is(t)の対数をガウス過程として直接モデル化することで、モデルを再パラメータ化します：
C(t)=logIs(t)-μs,�(�)=log��(�)-��、
2.10
ここで，C(t)は適切に選択されたゼロ平均ガウス過程であり，μsは感受性初期条件Is0に関係する感受性切片であり，次のようになる：
μs=logIs0-C(0).��=log��0-�(0).
2.11
この定式化は、主に切片パラメータとガウス過程との結合を緩め、サンプリングを高速化するために使用する。このことから、b(t)とlog Ir(t)を次のように計算することができる。
b(t)=ddtC(t)+γs�(�)=dd��(�)+�� です。
2.12
とし
logIr（t）＝C（t）＋μr＋∫t0γs dτ-∫t0u（τ）γT dτ-∫t0（1-u（τ））γU dτ＝C（t）＋μr＋∫t0γs-u（τ）（γT-γU）γU dτ＝C（t）＋μr＋（γs-γU）t-（γT-γU）∫t0u（τ）のdτ。 log��(�)=�(�)+�+∫0�� d�-∫0��(�)� d�-∫0�(1-�(�))�� d�=�(�)+�+∫0��-(��-��)-�� d=�(�)+�+（�-�）�- （��）0�（�）d�...�������(�)
2.13
もう一度、抵抗軌跡切片μrについて同じ推論を行い、Ir0と次のように関連付ける。
μr=logIr0-C(0).��=log��0-�(0) となります。
2.14
(d/dt)C(t)は、関連する共分散カーネルが十分に滑らかである限り、存在することに注意してください。微分可能な軌道を持つフルランクのガウス過程をメッシュ全体で評価することは，計算量がO(n3)（nは格子点数）となるため，法外なコストがかかる．このような高い計算コストは、このモデルを実現不可能なものにしてしまうでしょう。その代わりに、[34]で紹介されたフレームワークに基づいて、C(t)の低ランク表現に取り組みます。これは、C(t)の低ランク投影の表現につながり、C^(t)�^(�)で示される。
C^(t)=∑ j=1mSRBF( jπ2L---√;ρ,α)1L--√sin( jπ2L(t+L))f j�^(�)=∑ �=1�RBF( �2�;�,�)1�sin( �2�(�+�)) � �となる。
2.15
となり
ddtC^(t)=∑ j=1mSRBF( jπ2L---√;ρ,α)1L--√ jπ2Lcos( jπ2L(t+L))f j.dd��^(�)=∑ �=1�RBF( ��2�;゜, ゜)1� �2�cos( ���2�(゜+゜)) � �とする。
2.16
これにより、ガウス過程事前分布の評価複雑度は O(n3) から O(nm) に減少します。Lとmは先験的に指定する必要がある近似パラメータである（詳細は[34]を参照）。実際には、パラメータL = 6.5、m = 60のHilbert space Gaussian process (HSGP)近似を使用しました。これらの近似パラメータは、[34]に従って使用した長さスケールの事前分布の99%区間に適切である。ここで、fjは標準ガウス分布に従う独立かつ同一分布の確率変数、SRBF( - ; - , - )はRBFカーネルに対する適切なスペクトル密度、ρはカーネルの長さスケール、αはカーネルの限界標準偏差である[34]。
病原体動態モデルのパラメータをθθ＝（γs，γU，γT，Is0，Ir0，C＾（t））��＝（��，���，��0，�＾（�））で表す。ここで、モデル事後値 π(θ,α,ρ,f1:m∣gs,gr)�(�,�,�1:�∣��,���) を、適宜 t への依存性を抑えながら因数分解します：
π(θ,α,ρ,f1:m∣gs,gr)∝π(gs∣λs)π(gr∣λr)π(λs∣θ)π(λr∣θ)π(θ,α,ρ,f1:m) �(�,�,1:����,��)∝�（ �∣�）��（�� ∣∣∣�）� �(�� 1 :�) �（����） �（ ゜゜��）.
2.17
最初の 2 項は式(2.7)の合体尤度を用いて計算される。第3項は式(2.6)、(2.10)、(2.12)の組み合わせで与えられる。第4項は、式(2.6)、(2.12)、(2.13)を組み合わせることによって得られるものである。最後に、最後の項は次式で与えられる。
π(θ,α,ρ,f1:m)=π(C^(t)∣α,ρ,f1:m)π(γT∣γs)π(γU∣γs)π(γs)π(Is0)π(Ir0)π(α)π(ρ)π(f1:m),�(�,�,�,�1:�)=�(�^(�)∣�,�,�1:�)�(��∣��)�(��∣��)�(��)�(��0)�(��0)�(�)�(�)�(�1:�),
2.18
ここで、第1項はガウス過程（式(2.15)、(2.16)）で与えられ、残りの項は以下に示す事前分布に対応している。
2.5. 事前分布の選択とパラメータ化
このモデルは表1にまとめた事前分布でパラメータ化されている。データはγsの値についてあまり情報量が多くないと予想される。そのため、このパラメータにはかなり情報量の多い事前分布を設定し、先験的に知っていなければならない推測値γを中心とする。通常、95%区間で50%以上の相対変動を含むσ=0.3の値を使用します。σの値が大きいほど、幾何学が複雑になり、それに伴って事後分布のサンプリングも複雑になる。γTとγUはそれぞれ、耐性系統を対象とする焦点の抗生物質、または別の抗生物質で処理した場合の耐性系統の回復率を表す。γsを中心とし、正の値のみに切り捨てられた正規分布が自然な選択である。標準偏差を0.3γとすることで、重量の99%以上99%未満が2γ*以内に収まり、ありえないほど大きな変動が起こりにくくなる。このような大きな変動は、非常に急速な選択的掃討や絶滅につながるため、ここではほとんど関心がない。回復速度γTとγUは、式(2.5)を用いて、回復の絶対変化、つまりフィットネスパラメータに関連付けられる。γU＞γsは、耐性系統が感受性のある抗菌剤で処理された場合に回復が速くなることに対応し、したがって耐性のコストとなる。γT＜γsは、耐性系統が目的の抗菌薬で処理されたときに回復が遅くなることに対応し、したがって耐性の利点となる。代わりに、事後確率質量がγU < γsまたはγT > γsの割合が大きい場合、その結果は耐性コストまたは耐性ベネフィットが有意に存在しないことと一致すると結論づけられる。ρの事前分布は、質量の約1％がρ＜0.2の値に、約1％がρ＞2の値になるように選択した。下限はオーバーフィッティングを避けるため、上限はデータの範囲を超え、データから知ることができない長さスケールを抑制するために選択された。
表1.
モデルで使用したパラメータとプリオールの概要。
インラインで見るポップアップで見る
実際には、サンプリングアプローチを選択したため、γUとγTを制約のない空間でパラメータ化する必要があり、理想的にはγsへの依存性も弱める必要があります。そのために、パラメータq~U�~�とq~T�~�を導入し、γUとγTをこれらの決定論的変換と定義する：
andγU=log(1+exp{q~U+logγs})γT=log(1+exp{q~T+logγs}).⎫⎭⎬��=log(1+exp{�~�+log��})and�⎪�=log(1+exp{�~�+log��｝.
2.19
この変換に伴う尤度のヤコビアン調整は、次のように比例します。
|detJq|∝(1+exp{-q~T-logγs})-1(1+exp{-q~U-logγs})-1.|det��|∝(1+exp{-�~�-log�})-1(1+exp{-�~�-log�})-1.
2.20
2.6. 計算上の実装
式(2.17)の事後分布は高次元の分布であり、多くのパラメータが高度な相互依存性を持つことが予想されます。この分布からサンプリングするために、Stan [35]で利用できるHMCサンプラーである動的ハミルトニアンモンテカルロ(HMC)を使用します。HMCはマルコフ連鎖モンテカルロ法であり，エネルギー保存の特性を持つため，高い受入率を維持しながら，個々の状態間で大きなステップを踏むことができる．このため、微分可能な尤度を持つ中程度の高次元の事後分布からのサンプリングが効率的に行え、しかも反復回数が非常に少なくて済む。このモデルと推論方法はRパッケージとして実装されており、https://github.com/dhelekal/ResistPhy/。すべての結果は、ウォームアップに2000回、サンプリングに2000回の反復で、4つのチェーンを使用しました。すべてのモデルパラメータとすべての解析において、バルク有効サンプルサイズ（bulk-ESS）は常に500より大きく、すべてのRˆ�^統計は1.05より低く[36]、混合に問題がないことを示す値でした。また、少なくともサンプリング段階では発散遷移がないことを確認した。
2.7. シミュレーションと実データセットの使用
すべてのシミュレーションにおいて、式(2.1)の多系統SISの確率的な離散状態空間版を使用する。システムはタウリーピング[37]を用いてシミュレートされる。より具体的には、3つの系統を19年間かけてシミュレーションするシナリオを考える。2つの系統は感受性が高く、抗生物質の使用量の変動に影響されないように設定され、1つは耐性があるように設定される。最初の系統は、観測されていない集団の大部分を表すことを目的としており、したがって、はるかに高い有病率で開始されるように設定されています。2つの系統の軌跡を条件として、有効人口サイズNe(t)を変化させたキングマンの合体で系統樹をサンプリングし、軌跡を条件として式（2.6）に従います(28)。シミュレーションのパラメータは、耐性系統が102から104の間のオーダで有病率に達するように、妥当な範囲のもっともらしい行動を一貫して提供するように選択された。
CDC Gonococcal Isolate Surveillance Project（GISP）により、2000年から2013年の間に合計1102個のゲノムを収集した[38]。PhyML[39]を用いて最尤系統図を計算し、ClonalFrameML[40]を用いて組換えを補正し、BactDating[24]を用いて年代を算出した。この年代測定された系統樹は、以前、隠れた集団構造の解析で使用されたものと同じである[41]。1988年から2019年の間にGISPの参加者の間で淋病の治療に使用された主要な抗菌薬の分布は、https://www.cdc.gov/std/statistics/archive.htm で入手できるGISPレポートから得た。シプロフロキサシンとオフロキサシンの使用は、単一のフルオロキノロンカテゴリーに統合されたことに留意されたい。シミュレーションと実データセットの分析で使用したすべてのデータとコードは、https://github.com/dhelekal/ResistPhy/tree/main/run で入手可能です。結果
3.1. 単一のシミュレーションデータセットの詳細解析
このモデルの性能を検証するために、まず、図1に示すように、母集団サイズN（t）、感染率β（t）、抗菌薬使用関数u（t）が過去20年間変化する3系統の確率的SISによるシミュレーションを行う。最初の2つの系統は感受性が高く、抗菌薬使用量の変動に影響されないが、3番目の系統は耐性が高く、影響を受ける。最初の系統は感受性系統の大部分であるため、観察されないままになっています。残りの2つの系統は、それぞれ感受性のある系統と耐性のある系統で、観察された系統を表しています。感受性系統の1日あたりの回復率はγs = 1/60、qUに対する抵抗性の適合コスト = 1.25、qTに対する抵抗性の適合利益 = -2.7に設定された。これら2つの観察された系統のそれぞれから、200枚の葉を持つ年代別の系統樹がシミュレートされた。サンプリング日は最初の6年のうちの1年にランダムに割り当てられ、特定の年が選ばれる相対確率はその年の総有病率に比例した。このシミュレーションデータセットに対して推論を行った。その痕跡を電子補足資料の図S1に、カーネルパラメータの事後分布を電子補足資料の図S2に示す。感受性系統と耐性系統の有病率と再生産数R(t)を図2に示す。予想通り、推論値はシミュレーションで使用した正しい値を踏襲していた。図3に示すように、感受性系統の回復率γs、抵抗性のコストと利益qU、qTの推論値も、正しい値に近いことがわかった。γsの事後分布は、正しい値1/60を中心とする事前分布とほぼ同じであり、このパラメータについてデータが情報を持たないことを反映しており、情報量の多い事前分布を用いることの重要性が強調された。qUとqTの推定値の間には強い負の相関があった。この2つのパラメータは、抵抗性系統と感受性系統の全体的な適合度において正反対の役割を果たすので、予想通りである。しかし、qUとqTの推定値の範囲はそれぞれ1以上と1以下であり、平均値を1とする対数正規分布に反していることから、抵抗性に関連するコストと利益の両方を検出することができた（図3）。最後に、時間A(t)を経た祖先系統の数の事後予測分布[42]を計算し、入力した系統データと比較した（電子補足資料、図S3）。データと事後予測分布は類似しており、シミュレーションと推論に同じモデルを用いていることから、予想されるように、モデルとデータとの適合性が高いことが示された。
キャプション
図をダウンロードする
新しいタブで開く
パワーポイントのダウンロード
キャプション
図をダウンロードする
新しいタブで開く
パワーポイントのダウンロード
キャプション
図をダウンロードする
新しいタブで開く
パワーポイントのダウンロード
3.2. 複数のシミュレーションデータセットによるベンチマーク
フィットネスコストと抵抗の利益の値を変化させた以外は、上記、図1に示したのと同じ条件でシミュレートしたデータに対して、推論手法を同様に適用することを繰り返した。フィットネスコストqUは1から1.2まで直線的に増加し、フィットネスベネフィットqTは1から0.5まで直線的に減少する、合計50個のシミュレーションデータセットが作成され、分析されました。これらのシミュレーションにおける感受性系統と耐性系統の有病率は、電子補足資料の図S4に示されている。推論の結果は図4に示されており、ほぼすべてのケースで、事後の95%信頼区間がシミュレーションで使用した適性コストと抵抗性の利益の正しい値をカバーしていることが示されている。
キャプション
図のダウンロード
新しいタブで開く
パワーポイントのダウンロード
3.3. 米国におけるフルオロキノロン耐性N.gonorrhoeaeへの適用
淋菌におけるフルオロキノロン耐性のコストと利益を推定することで、我々のモデルと推論フレームワークの使用を実証します。CDC GISP [38]が2000年から2013年の間に収集した1102個のゲノムに基づいて、ClonalFrameML [40]を用いて組換え補正されたツリーを構築し、BactDating [24]を用いて日付を決定しました。この系統樹には2つの主要なフルオロキノロン耐性系統が存在するため[38]、比較研究を行うことにした。2つのフルオロキノロン耐性系統と1つのフルオロキノロン感受性系統は、他の関連する抗生物質に対する耐性プロファイルの類似性に基づいて選択されました。抗生物質の使用データと3つの系統の耐性プロファイル（図5）を見てみると、1995年以降に一次治療としてかなりのレベルで使用されていた抗菌薬の耐性プロファイルが一致していることがわかります。そのため、この年を解析開始日（tmin = 1995）、最後のゲノムを収集した日を終了日（tmax = 2013）としました。なお、感受性系統のうち、セフィキシムに対するde novo耐性獲得を示したサブクレードは削除しています。感受性系統の1日あたりの回復率の事前平均は、過去の淋病のモデル化研究[10,43,44]に基づいてγ* = 1/90 に設定された。
キャプション
図をダウンロードする
新しいタブで開く
パワーポイントのダウンロード
このデータセットに対して推論を行った。トレースは電子補足資料の図S5に、カーネルパラメータの事後分布は電子補足資料の図S6に示されている。図6は、2つの耐性系統の事後潜伏感染ダイナミクスの要約を示し、電子補足資料の図S7は、感受性系統の同じものを示す。2つの耐性系統は同様のダイナミクスを示し、フルオロキノロンの使用が減少した瞬間に相当する2007年頃に有病率がピークに達している（図5）。図7は、両耐性系統の耐性パラメータqUとqTの周辺分布と結合事後分布を示したものである。これは、qTとqUがそれぞれ高い事後確率で1以下と1以上に局在することから、両系統ともフルオロキノロン耐性にコストとベネフィットの両方が存在することと矛盾しない。このことは、両系統が同じような選択圧に直面し、フルオロキノロン耐性に関連する体力コストを克服するために適応に成功したとは思えないことを示しています。モデルがデータを適切に説明できることを確認するために、事後予測アプローチを使用した[42]。時間経過に伴う祖先系統の機能A(t)の事後予測軌跡をシミュレーションしたところ、系統データが示唆するものと非常に似ていることがわかった（電子補足資料、図S8）。
キャプション
図をダウンロードする
新しいタブで開く
パワーポイントのダウンロード
キャプション
図をダウンロードする
新しいタブで開く
パワーポイントのダウンロード
系統間の完全競争の仮定において、抵抗性系統が確立できず、その割合が十分に速く減少することを保証したい場合、減衰係数c > 0を固定し、抵抗性系統の成長率が感受性系統の成長率よりc単位低くなるようにする、すなわちrs（t）- rr（t）>cとなるようにします。成長率は回復率の誤仕様の影響を受けにくいので、成長率を使うことにした。系統が同じ伝達率関数b(t)を持つ場合、この条件はγs(t) - γr(t)>cと等価であり、式(2.3)のγr(t)の定義を使えば、u(t)qT + (1 - u(t))qU >cと等価である。図8に示すように、感受性系統と両耐性系統のそれぞれの成長率の差はcを超える。耐性系統が感受性系統よりも低い適合度を維持することを95％確実にするためには、耐性系統1、2について、それぞれ感染者の約20％20％、15％15％以上にフルオロキノロンを処方してはいけない。
キャプション
図をダウンロードする
新しいタブで開く
パワーポイントのダウンロード考察
ある抗生物質に耐性を持つ細菌病原体の系統は、同様の感受性を持つ系統と比較して、フィットネスコストとフィットネスベネフィットの両方が発生する [8] 。抗生物質が広範囲に使用されている場合、その利益はコストよりも大きくなる可能性が高い。この場合、耐性系統は感受性系統に対して選択的に有利であり、したがって、より速い速度で成長する。逆に、抗生物質がほとんど使われない、あるいは全く使われない場合、利益はコストよりも小さくなり、耐性系統の頻度が減少する可能性があります。したがって、これらのパラメータを推定することは、耐性菌の増加を招くことなく抗生物質を処方する方法を決定する上で、最も重要なことである[9]。ここでは、流行ダイナミクスと系統力学の関連性を示した過去の研究 [13,19,20,28] に続き、ゲノム配列データと抗生物質の処方に関するデータを組み合わせて、この目的に利用できることを示した。感受性系統と耐性系統の系統力学的軌跡を比較し、それらを抗生物質使用量の既知の関数と関連付けることで、耐性の適合コストと利益に対応するパラメータを個別に推定できることを示す。特に、N. gonorrhoeaeゲノムの大規模な公開コレクション[38]を再解析した。その結果、フルオロキノロン系抗菌薬に耐性を持つ淋菌の2つの系統について、これらのパラメータを推定することができ、両者で同様のコストと利益を推定することができました（図7）。この知見をもとに、フルオロキノロン系抗菌薬のスチュワードシップに関する提言を行うことができました（図8）。
我々の手法のインプットとして、感受性系統と耐性系統の両方について、日付の入った系統樹が必要である。BEAST[22]やBEAST2[23]などの配列アライメントから、あるいはtreedater[45]やBactDating[24]などの未作成の系統樹から、いくつかのソフトウェアツールを使ってこれを作成することができる。このような日付のある系統樹を構築するには、サンプリング期間中に集団が測定可能な進化を遂げているか [46,47]、または分子時計速度の過去の推定値 [48] が必要です。本手法が必要とするもう一つの入力は、関連する時間枠と地理的な場所での抗生物質使用関数である。これは、すべての歴史的背景で常に利用できるとは限りませんが、これらのデータを収集するための努力がますますなされています[49]。最後に、本手法では、感受性系統の回復率に関する情報量の多い事前分布を必要とする（表1）。これは、多くの類似のコンパートメント型流行モデル[50]と同様に、データから特定できないのが普通だからである。この事前分布は、研究対象の感染症によって、また既存の科学文献に基づいて慎重に選択する必要がある。
我々の推論方法は、定義が明確で比較的単純な疫学モデル（式(2.4)）に基づいており、これは多くの仮定を立てることを意味し、その妥当性は分析を行う前に検討されました。このモデルでは、多系統の病原体の動態が、十分に混合された宿主集団における人から人への感染によって引き起こされ、有意な集団構造が存在しないため、系統間の競争が完全であると仮定しています。また、個体は感染するとすぐに感染力を持ち、回復するまでその感染力は一定であり、回復後は免疫力を得ることなく再び感染力を持つようになると仮定している。このように比較的強い仮定が並んでいるため、実際の感染症には適用できないように思われるかもしれないが、推論が可能なモデルを得るためには必要である。さらに、これらの仮定に違反しても、推論の結果が無効になるとは限らない。例えば、感染が免疫を引き起こす場合、感受性個体の数S(t)を効果的に減少させるが（式(2.2)）、本モデルではこの数は一定であるとは仮定していない。実際、宿主集団の大きさN(t)と感受性個体の数S(t)は、関数b(t)のパラメータ化の一部として統合されているので（式(2.4)参照）、付与された免疫が研究対象の全系統に当てはまる限り、推論は頑強である。同様に、非構造化集団の仮定は、アメリカ全土の淋菌への適用を含め、問題があると思われるかもしれませんが、小規模な局所的なアウトブレイク以外では、解析に利用できるゲノムは感染集団全体から疎にサンプリングされます[51]。このような状況では、実際の感染者数ではなく有効な感染者数を考慮する限り、宿主の集団構造が系統力学に及ぼす影響は軽微であると考えられる[52,53]。
解析中の系統データと我々のモデルの適合性は、事後予測分布チェックを用いて検証できる（電子補足資料、図S3およびS8）。これらのテストが失敗した場合、またはモデルの仮定が不適切であると考えられる場合、解決策として、抵抗性のパラメータを直接推論しない代償として、年代物の系統を後処理する他の方法 [25] に頼ることができるが、仮定は少なくなる。代替アプローチとしては、異なる系統における分岐パターンの違いを検出するノンパラメトリック手法[41,54]や、根本的な疫学的要因ではなく、病原体の集団サイズの成長という観点からパラメータ化した手法[15,55]がある。しかし、我々のモデルベースのアプローチは一般的で柔軟性があるため、https://github.com/dhelekal/ResistPhy/ で公開されている我々のソフトウェアの実装を使用すれば、多くの環境で適用できることが期待されます。この手法を多くの細菌病原体に関するますます大規模になるゲノムデータベースに適用することで、抗生物質の使用量と耐性との間の正確な関連性を定量化し、将来の抗生物質処方戦略の設計に必要なエビデンスベースを提供することができると考えています[9,56,57]。
データ入手の容易性
本研究で使用したすべてのデータとコードは、GitHubのデジタルリポジトリ（https://github.com/dhelekal/ResistPhy/）から入手可能です。
データは電子補足資料[58]で提供されています。
著者の貢献
D.H.：概念化、調査、方法論、ソフトウェア、執筆-原案、執筆-レビュー、編集、M.K.：概念化、資金獲得、調査、監督、執筆-レビュー、編集、Y.H.G.：概念化、調査、執筆-レビュー、編集、X.D.：概念化、資金獲得、調査、方法論、監督、執筆-原案、執筆-レビュー、編集．
すべての著者が出版を最終的に承認し、その中で行われた仕事について責任を負うことに同意している。
利益相反の宣言
我々は、競合する利害関係がないことを宣言する。
資金提供
National Institute for Health Research (NIHR) Health Protection Research Unit in Genomics and Enabling Dataからの資金提供を謝意を表します。この研究は、英国工学・物理科学研究評議会（EPSRC）の助成金No. EP/S022244/1（EPSRC Centre for Doctoral Training in Mathematics for Real-World Systems II）の助成を受けました。
脚注
電子補足資料は、https://doi.org/10.6084/m9.figshare.c.6673522 でオンライン公開されています。
© 2023 The Authors.
原著者および出典を明記することを条件に、無制限の使用を許可するクリエイティブ・コモンズ表示ライセンス http://creativecommons.org/licenses/by/4.0/ の条件の下で、王立協会が発行したものです。
全文を見る pdfをダウンロードする
戻る
インターフェイス
本誌について
著者特典
投稿
お問い合わせ先
トランスフォーメーショナルジャーナル
購入情報
ジャーナルメトリクス
図書館への推薦
検索ヘルプ
ロイヤルソサエティ出版
ジャーナル
歴史的背景
オープンアクセス
オープンサイエンス
出版ポリシー
利用許諾
コンファレンス
ビデオ
ブログ
アカウント管理
ご利用条件
個人情報保護方針
クッキー
ロイヤルソサエティ
ロイヤルソサエティについて
お問い合わせ
フェロー
イベント情報
助成金・スキーム・賞
トピックス・政策
コレクション
会場レンタル
トップページへ戻る
著作権 © 2023 英国王立協会
この記事が気に入ったらサポートをしてみませんか？