皆さん、こんにちは。ベーリンガーインゲルハイムで主席統計学者として働いているマルコ・クンツェルマンです。JMPのクリス・ゴットウォルト氏とは、DOEとクロマトグラフィーデータを組み合わせるというテーマで素晴らしいコラボレーションをしています。まだ作業は進行中ですが、このプロジェクトは既に最新の結果をお見せできる段階に達していると考えています。なお、プレゼンテーションで紹介されている多くの手順は依然として手動で行われていますが、JMPの将来のバージョンでは自動化される予定です。
このプロジェクトの背後にある動機は、言うまでもなく、効率向上とコスト削減が産業界における基本的な目標であるということです。また、近年、規制当局や規制ガイドラインは、多変量設計空間（DSP）の特性評価と調査というテーマをますます推進しています。その一例として、ICH Q8と最近発表されたICH Q14が挙げられます。これらの規制当局の要件を満たす確立された方法の一つが、実験計画法（DOE）です。
2024年6月に発効したICH Q14の発効により、DOEを用いた統計モデリングが分析開発に活用されるようになりました。しかし、クロマトグラムは通常非常に複雑であり、分析法開発に統合しようとする人々は、クロマトグラムから特徴を抽出し、UE入力パラメータの関数としてモデル化するのが一般的です。
どういう意味でしょうか？特徴抽出とは、クロマトグラムから得られるピーク高さやピーク保持時間、あるいは2つまたは3つのピークが互いにどれだけ分離されているかを表す分析分解能といった情報です。特徴抽出の問題点は、そもそもピーク検出のためのソフトウェアが必要なため、特徴抽出が面倒な作業になることです。そして通常はピーク積分が必要で、その後、分析ソフトウェアからJMPにこれらの情報をすべてエクスポートする必要があります。JMPでは、後から通常のDOE法を用いてモデリングを行うことができます。
さらに、抽出される特徴それぞれに1つのモデルが必要です。先ほどの例では、9つの異なるピークが見られましたが、ピークの保持時間を記述するだけで9つのモデルが必要でした。さらに、分解能、ピーク高さなどを記述するためにも9つのモデルが必要になります。これらすべての個々のモデルを最適化しようとすると、最終的に解釈が非常に困難になることは想像に難くありません。なぜなら、保持時間、ピーク高さ、分解能などの情報しか得られず、全体像を把握できないからです。完全なクロマトグラムを見ることはできないのです。
これが、このプロジェクトで私たちが達成したかった基本的な目標です。JMPの機能データエクスプローラー内で、すぐにピーク検出ができるようにしたいと考えました。そして、もちろんJMP内で、完全なクロマトグラムの数値パラメータ化を行いたいと考えました。すべてのピークを含む完全なクロマトグラムに対して、直接モデリングを適用する必要があります。そして最後に、もちろんJMPのインタラクティブ機能を活用したいと考えました。つまり、プロファイラー内ですぐにインタラクティブなクロマトグラムを作成し、異なるDOE係数設定がクロマトグラム自体にどのような影響を与えるかを実際に確認できるようにしたいと考えました。
クリスと私は、この1年間に扱ってきた2つの異なるデータセットの例を示すことで、私たちのプロセスを最も分かりやすく説明できると考えました。まずは、ガスクロマトグラフデータを用いた最初の例から始めましょう。ガスクロマトグラフデータ（略してGC）は、広く用いられている分析手法です。簡単に言うと、サンプルを蒸発させると、ガス相において、移動相とキャピラリー自体にある固定相との間で相互作用が起こります。これらの相互作用効果は、分析サンプル中の成分ごとに異なります。これに基づいて、これらの異なる成分を互いに分離し、最終的に検出器で検出することができます。
この例では、超高速GCと呼ばれる別のGC装置を使用しました。分離原理は通常のGCと同じです。しかし、ここでの主な違いは、加熱速度が大幅に向上していることです。これにより、完了時間が1分以下でクロマトグラムを作成できるようになりました。これはもちろん、このプロジェクトにとって非常に有益です。多くの異なる設定を簡単に実行でき、時間をあまり無駄にすることなく実行できるからです。
ここでのDOE実験は、4つの異なるDOE因子に基づいています。つまり、開始温度、終了温度、保持時間、そしてデルタタイムです。こうした情報をデータテーブルだけでなく、右側の温度プロファイルのフォームにも表示しています。これが初期温度です。この初期温度を変化させ、保持時間を変えながら温度を維持しました。
異なるΔtを適用することで、異なる加熱速度と、それと同時に異なる最終温度（T₂）を適用することができました。この超高速GCでは多くの実験を容易に実行できるため、システムから得られるすべての情報を実際に得るために、81回の実験と7回の中心点実験を含む完全実施要因計画を採用することにしました。
そのため、DOEでは、4つのDOE因子がクロマトグラムの9つの異なるピークに与える影響を調べました。下の部分はGCクロマトグラムの一例で、9つの異なるピークの成分が生データとして表示されています。最終的に、88回のDOE実験すべてについて、このクロマトグラムが得られると考えられます。
ここからはクリスにバトンタッチします。クリスは、最終的に合成クロマトグラムを作成するために、クロマトグラムからどのように情報を抽出したかを詳しく説明します。クリス、お願いします。
ありがとう、マルコ。これは、プロファイラーで合成または予測されたクロマトグラムを4つのDOE因子とともに視覚的に表現した究極の例です。合成クロマトグラムは左側に黄色の背景で表示されています。モデル化する必要があった特定のピークが9つあります。
これからご紹介するアプローチでは、プロファイラーコントロールを用いた分かりやすい視覚化を作成し、分析者は最適なクロマトグラムが得られるDOE係数設定を調整できます。今回のケースでは、ピークが最も分離している、またはピーク分離度が高いクロマトグラムです。設計空間の中央にあるこの設定（複数のピークが密集している状態）から、本研究で得られた最良の分離度を達成したこれらの設定まで、様々な設定が可能です。
これを実現するために、克服すべき課題が数多くありました。まず、不要な溶媒ピークとブラインドベースラインランを除去する必要がありました。そこから、JMP Pro 18の機能データエクスプローラーのピークモデリング機能を用いてピークを識別し、ラベル付けを行いました。すべてのクロマトグラムについて、各ピークのピーク高、半値幅、溶出時間を抽出しました。次に、これらをGenregを用いてモデル化し、得られた予測式を全体の合成クロマトグラム式にまとめました。
昨年この取り組みを始めた当初は、多くの手作業が必要でした。Ryan ParkerとClay Barkerは、JMP Pro 18でこのプロセスを効率化するために多大な努力を払いました。デモをご覧いただければわかるように、まだ多くのステップがありますが、数日かかっていた手作業が約1時間に短縮されました。今後のJMP Proリリースでは、さらにプロセスを効率化していく予定です。
まず、いつもの方法でFDEに入ります。その前に、各バッチの元の開始時間と終了時間を計算するために、tabulate関数を使う必要がありました。これは、データのアラインメントを解除し、ブラインドラン用のモデルをスケーリングして、各クロマトグラムから差し引くために必要です。ここには生のクロマトグラムが表示されています。画像の大部分を占める溶媒ピークを除去したいので、Align Maximum関数を使用します。
これで、溶媒ピークの後、かつ最初の注目ピークの前の、アラインメントされた時間を見つけることができました。十字線ツールを使ってこれを行うと、0.26がすべての溶媒ピークと最初の注目ピークの間の時間であることがわかります。0.26より前の観測値はすべて除外します。これで、データから注目ピークだけが得られました。これを新しいデータセットとして保存し、作業内容を明確にします。
ベースラインのブラインドランは、製品を添加しない状態で実施されました。このシグナルはモデルに悪影響を与え、結果にバイアスを与えるため、データから除去する必要があります。ここでは、ブラインドランのみをpスプラインでフィッティングし、サマリーを保存しました。次に、溶媒ピークを除去したデータに式をコピーしました。次に、ブラインドランの式を使用してクロマトグラムを手動で整列させ、絶対時間スケールに戻したブラインドランの式列の強度を作成する際に、ブラインドランオフを差し引きました。
クロマトグラムデータが目的の形式で取得できたので、FDEに戻りましょう。JMP Pro 18の新しいピーク検出機能を使用します。まずは自動ピーク検出を実行します。これにより、データ内の多数のピークが検出され、検出された候補ピークごとに範囲と中心位置が示されます。「サマリー」タブをクリックすると、最も重要な導出ピーク指標を確認できます。
ここでは、AUC、分離能、ピーク高さ、その他のピークサマリー指標を確認できます。多くのクロマトグラムのサマリーに9つ以上のピークが含まれているため、問題があることに気付きました。しかし、実際にはスペクトルごとにピークは9つしかありません。
そこで、ピーク検出ユーザーインターフェースに戻り、9つ以上のピークを持つスペクトルを調べます。少し調べてみると、高さが500未満のピークと思われるものはすべて、データ内の小さな偶然によるものであることがわかります。ここでは、こうした偽ピークの一例をハイライト表示しています。自動ピーク除去ツールを使ってこれらのピークを除去するには、ツールをクリックし、有効なピーク高さの下限値として500を入力します。
500未満のピークを除外することで、合成クロマティックRAMモデルの作成に必要なピーク要約統計量のクリーンかつ正確な抽出が可能になりました。これらの要約を新しいテーブルに保存し、一般化回帰分析を用いて、抽出したすべての指標に対数正規分布の2次応答曲面モデルを当てはめます。すべての指標が正の値であるため、対数正規分布を使用します。
実際には、他のモデリング演習と同様に、モデルが適切に適合しているかどうかを確認する必要があります。ここでは省略しますが、満足のいく結果が得られたら、コントロールボタンを押したまま予測式を保存するためのオプションを選択することで、すべてのレスポンスの予測式を一挙に保存できます。これで、DOE係数を用いてクロマトグラムの各特性を予測するモデルができました。そして、各ピークの位置、半値幅、幅、高さについて、個々のDOEモデルを組み合わせた式列を用いて、合成クロマトグラムを計算する予測式を作成しました。
モデルを検証するために、赤で描かれた合成クロマトグラムと青で描かれた実際のクロマトグラムの重ね合わせプロットを見てみましょう。全体的に見て、モデルは完璧ではありませんが、実際のスペクトルの主要な特徴をかなりうまく捉えているように見えます。ここからは、プロファイラーを使用して、計画の中心から、すべてのピークに対して製造プロセスの成分を可能な限り適切に分離する分析手法を与える設定まで、ピークの全体的な分解能を向上させる因子設定を見つけます。
これで、分析法開発に携わる方々に分かりやすくプロセスを設定する方法を示した、シンプルで分かりやすいビジュアライゼーションが完成しました。それでは、Marcoにバトンタッチしましょう。
クリス、本当にありがとうございます。クリスは、インシリカクロマトグラムがGCデータとどのように連携するかを示してくれました。この手法は、キャピラリーゲル電気泳動という別の例にも既に適用しています。CGEは、簡単に言えば、分析機器の一種です。ただし、分離技術が少し異なり、ピーク形状も異なりますが、これについては後ほど説明します。CGE法開発の主な動機は、2つの異なるピークを分離する必要があったことです。率直に言うと、ここではピーク1とピーク2を分離しています。
ピーク1はピーク2に比べて非常に小さく、しかもこれら2つのピークは互いに大きく重なり合っているため、これは非常に難しいことがわかります。ここでも、これら2つの成分を分離するための最適な設定を見つけようと試みます。今回は、最終濃度を除いて主に3つの異なるレベルで、5つの異なるUE係数を設定しました。最終濃度は4つの異なるレベルで調査しました。
通常、ここで全ての組み合わせを調査する場合、完全実施要因設計では324回の実験が必要になります。これが、最終的にI最適設計を採用した主な理由です。これにより、わずか64回の実験設計を作成したため、約80%の実験回数を削減できました。この設計は、主効果、2因子交互作用、二次効果を解析できるだけでなく、3因子、4因子、5因子の交互作用効果、そして主効果と二次効果を正確に組み合わせて解析することができました。
なぜそうしたのでしょうか？GCの例では、多くの分析DOE要因が相互作用していることを実際に調査しました。そして、これをさらに詳しく調査し、最終的に高い精度と正確性を得るためには、モデルがどれほど複雑になる必要があるかをより深く理解したかったのです。
すでに述べたように、通常はクロマトグラム全体をモデル化し、特徴量をまず抽出することも可能です。しかし、今回のケースではそれが全く不可能でした。2つのピークが非常に重なり合っているため、解析ソフトウェアはこれらの2つの成分の解析分解能を計算できなくなっていたのです。これが、クロマトグラムを数学的に記述する別の方法を見つける必要があったもう一つの理由です。
ここで何をしたのでしょうか？クリスが既に説明したように、ピークは当然ながら数学関数の形で記述できます。ピークの高さ、保持時間、ピークの幅などです。対称ピークは、この左側の式の形で記述できます。
残念ながら、CGEでは状況はより複雑になります。なぜなら、CGEは常に非対称なピークをもたらすからです。ここでも、2番目の成分、つまり2番目のピークは全く対称的ではないことがわかります。このピークの非対称性を記述できる別の式を見つける必要があります。
ここで文献で見つけたのが、この便利な式です。パラメータaは、ピークの非対称性を記述できるパラメータです。左側と右側を見ると、値が0でもピークが非対称であることがわかります。aの値が大きくなるほど、ピークの非対称性も増します。その他のパラメータは、基本的にピークの保持時間やピークの高さなどです。これは基本的に前の式と同じですが、少し複雑に見えます。
この式をJMPに転送するにはどうすればよいでしょうか？JMP内に新しい式列を作成しました。「列」→「式」と進み、「式エディタ」を開きます。左下の「パラメータ」→「新規パラメータ」をクリックし、必要なパラメータをすべて追加します。
そうすることで、緑の四角で囲んだ部分のように、2つのピークの式を作成できます。ピーク1については対称ピークを仮定した式、2つ目のピークについては、ピークの非対称部分も記述したいため、より複雑な式を作成します。また、ここには、2つ目のピークのh1、w1、t1、m2、d2、a2、H2といった、入力したすべてのパラメータも確認できます。次のステップでは、これらすべてのパラメータに適切な値を見つける必要があります。
そのためには、まず分析機器からすべての生データをJMPに読み込みます。ここでは、移動時間を表す列と強度を表す列が1つずつあります。個々のエレクトロフェログラムは、3つ目の列「Run」で区切られており、これはDOEのRun番号です。
その後、非線形プラットフォームの「解析」→「特殊モデリング」に進みます。ここでは、Y応答と強度を入力し、予測式には新しく作成した式列を使用します。そして、すべてを列ごとに分けます。ここに示したのは「実行」列で、個々のDOE実験がすべて含まれています。
その後、「OK」をクリックすると、ご覧のようなビューが表示されます。これは2回目のDOE実験で、少しズームインして26分から28分までのデータのみを表示しています。点は実際の生データです。実線はパラメータ化された関数で、パラメータ推定値は下図のようになっています。初期の値は当然ながら最適な値ではありません。
次のステップでは、JMPを実行して、生データに可能な限り適合させるために、個々のパラメータすべてについて適切な推定値を見つける必要があります。これは反復的なステップです。1回の反復を実行するたびに、最適な結果に徐々に近づいていきます。ここでは、この小さな図と4～5回の反復でこれを説明しようとしています。もちろん、実際には250回の反復があり、64個の実験すべてについてこのプロセスを繰り返す必要があります。
最適化ツールがすべてのパラメータについて本当に適切な推定値を見つけたかどうかを視覚的に確認することは非常に重要です。これが可能になるには、適切な初期値を見つける必要があります。そうでなければ、局所最適解に陥るか、アルゴリズム最適化ツールが全く収束しないという問題が生じます。もちろん、これはここでの問題です。適切な初期値を見つけるためには、アルゴリズムに下限と上限を与え、まずクロマトグラムを調べてピーク高さや保持時間などについて、少なくとも大まかに適切な推定値を得ることが理にかなっています。これは現時点では非常に面倒で手間のかかるステップですが、JMP 19では自動化される予定です。
すでに大変な作業でしたが、64回の実験のうち9回が失敗の淵に陥ったため、課題はさらに大きくなりました。つまり、2つのピークの移動時間は、いくつかのDOE因子の組み合わせにおいて予想外に長くなりました。この手法は60分後に記録とともに自動的に停止するため、9回の実験ではこれらのピークに関するデータは得られませんでした。これで2つの選択肢が生まれました。これらの実験を単純に失敗とするか、JMPに既に実装されている手法を用いて、このピークから得られる情報の少なくとも一部を活用するかです。
この手法はデータ打ち切りと呼ばれます。例えば、ピークの移動時間については、ピークは少なくとも60分後に現れることが分かっており、専門家は少なくとも60分から62分の間であるはずだと答えることができるので、この情報を利用できます。同じことは、ピークの高さや幅などの他のピークパラメータにも適用できます。
JMPでこの情報をどのように活用したのでしょうか？すべてのパラメータに対して、lowerとupperという新しい列を作成するだけです。この例では、最も簡単なパラメータM2、つまり2番目のピークの保持時間または移動時間を使用しました。ここではM2_lowerとM2_upperです。
データが得られたすべての実験において、下限成分と上限成分に同じ保持時間を入力すれば済みます。これは2回目、3回目、そしてそれ以降の実験でも同様です。ただし、移動時間が得られなかったピーク、つまり失敗限界値に達した実験では、そのピークの下限予測移動時間を使用します。実験1では60分、上限予測移動時間は62分です。a2、d2などの他のピークパラメータについても同様に処理を進めます。
次のステップは非常に簡単です。「モデルのあてはめ」プラットフォームの「一般化回帰」に進みます。M2_lowerとM2_upperを挿入し、もちろん、DOE因子、因子の組み合わせ、交互作用効果などをすべて入力します。「実行」をクリックするだけで、JMPは区間打ち切りデータの使用を自動的に認識し、それが正しいかどうかを確認します。ここでは「はい」をクリックします。これにより、すべてのデータに上限と下限の情報を使用できます。これは、単にデータを削除するよりもはるかに優れています。なぜなら、失われたデータに関する情報が得られており、それをモデリングに活用できるからです。
ピーク高さa2、非対称パラメータa2など、他のすべてのパラメータについても同様の手順で進めます。また、最初のピークについても同様の手順で進めます。次のステップでは、モデルのフィッティングとモデル選択（ここでは一般化最小二乗法プラットフォームのPruned Forward Selection法を使用）の後、「予測式を保存」を使用して、各パラメータの列を保存します。
次のステップは実は非常に簡単です。このプロセスをh1、t1、w1など、すべてのパラメータについて繰り返します。次に、パラメータを対応する予測式に置き換えるだけで、クロマトグラムの式を再構成します。
覚えていると思いますが、最初にピーク1とピーク2について複雑な式を使っていました。そして今、個々のパラメータを予測式に置き換えています。これで予測式にh1はなくなり、式列に入力されました。代わりにh1_upper予測式が表示され、この処理が他のすべてのパラメータに対して繰り返されます。こうすることで非常に長い式が生成されるため、画面上では収まりきらないように見えますが、最終的に完全な式がどのように続くかは想像できると思います。
これが結果です。JMPプロファイルでインタラクティブな方法で結果を確認してみましょう。次のステップでは、グラフプロファイルを開き、強度の式（ここで示したように再構成された式）を選択します。「中間式を展開」をクリックし、強度の式とY予測式を入力して「OK」をクリックします。
JMPで結果がどのように見えるか、簡単なデモをお見せしたいと思います。こちらはJMPのプロファイルプラットフォームです。左側には、緩衝液pH、DOE係数、緩衝液SDS、DOE係数など、実際に調査したすべてのDOE係数が表示されています。
興味深いのは右側の電気泳動図です。X軸は泳動時間、Y軸は強度を示しています。既に少し拡大表示していますが、これが主に35～47分の範囲のグラフを表示している理由です。しかし、2つのピークが互いにどれだけよく分離されているかを確認すると、さらに良く分かります。プロファイラーやDOE係数を直接変更し、クロマトグラム自体で何が起こっているかを確認できるようになりました。最終濃度を変えると、当然のことながらピークの高さも変化します。
しかし、分離電圧を変えると、これら2つのピークの移動時間も変化します。このプロセスをキャピラリー温度と緩衝液SDSで繰り返すことができます。よく見ると、2つのピークの分離が改善されている場合もありますが、この画像のように、実際には分離が不可能なほど悪化している場合もあります。しかし、緩衝液のpHを下げると、これら2つの成分の分離が大幅に改善されます。
これはほんの一部です。これでJMPの双方向性と合成クロマトグラムが組み合わされました。しかし、さらに一歩進めてみたいと思いました。ここで選択した設定でエレクトロフェログラムを少し拡大し、予測値と信頼区間の式を調整しました。信頼区間と予測区間の式をピーク1とピーク2の式に置き換えました。信頼区間の式には、ピーク1の成分だけが表示されるようになりました。
もう少し大きくしてみましょう。予測区間の式は、2番目のピーク成分のみを示しています。これで非常に興味深い結果が得られました。2つのピークが重なり合っているのが分かります。個々のピークがどこに現れ、どの程度離れているかも確認できます。
これで、ピーク1の成分が2番目のピーク成分のどこにあるかを正確に把握できるようになりました。これはもちろん、分離がどの程度適切に行われているか、そして最終的にどの設定が最適かを理解するのに大いに役立ちます。
モデリング後に必ず行うべきことは、実際の生データとモデルの予測値を比較することです。ここでは、予測値と測定値のプロットでこれを行いました。赤で示されているのが実際のデータ、青で示されているのがモデルの予測値です。ここでは、モデルのパフォーマンスがいかに優れているかを示すために、いくつかの実験を選択しました。
混乱しないでください。まずはRun 2から始めます。Run 1は失敗の瀬戸際だった実験の一つで、そのためピーク成分のデータがなく、モデル予測との比較もできていません。ここではあまり役に立ちません。しかし、Run 2からRun 11までの他の全ての実験では、モデルが実際にはかなり良好に機能していることがわかります。保持時間に若干の変化が見られる場合もありますが、ピーク高さは通常かなりよく適合しています。また、ピーク幅もかなり良好です。
ここではベースラインに少し問題がありますが、クリスと彼のチームは現在、このベースラインモデリング技術を用いてこの問題に真剣に取り組んでおり、改善に取り組んでいます。全体として、モデルのパフォーマンスを見ると、実データとモデルの間で非常に良好な予測が得られていると言えるでしょう。
CGEの例で最後に行いたいのは、もちろん、ピークを可能な限り分離するために最適なDOE係数の設定を見つけることです。通常は、2つのピーク成分の分析分解能を単純に計算します。保持時間、つまりこの場合は2番目のピークの移動時間から1番目のピークの保持時間を差し引き、それを1番目のピークと2番目のピークのピーク半値幅で割ります。
これらのパラメータのモデルはすべて、DOE係数の関数として既にモデル化されています。これらの情報をプラグインアンドプレイするだけで、調査したDOE係数の関数として分解能を持つ、いわゆるヘルパーモデルを作成できます。クロマトグラムの右上に、分解能が表示されています。その上の2行目には、DOE係数の関数として分解能が表示されています。
ここでJMPの機能を使って設定を最適化します。つまり、分解能を最大に最適化するように指示するだけです。境界条件として、専門家は移動時間が50分を超えないようにしたいと考えていました。そこで、M2予測式（既にある）を上限の50分に設定しました。次に、アルゴリズムを実行して全てを最適化し、結果を取得します。これがここに表示されているものです。これらのDOE係数設定により、2つのピーク成分を適切に分離でき、CGE法開発の問題を解決できました。
これは素晴らしい成果でした。まとめると、重要な点は、JMP内でピーク検出が可能になったことです。将来的には、クロマトグラムの数値パラメータ最適化もJMP内で直接自動化する予定です。現在はまだ手作業で、ご自身で行うことも可能ですが、それでも面倒な作業です。しかし、JMP 19ではこの点が改善される予定です。
すべてのピークを含む完全なクロマトグラムを直接モデル化し、その結果をJMPプロファイルで直接確認・視覚化できるという画期的な機能が追加されました。私の経験から言うと、これは各分野の専門家にとって、DOEの仕組みや手法の最適化方法を理解する上で大きなメリットとなりました。これはまさに、統計学者と研究室の現場、そして分析化学者との間の溝を埋める最先端技術と言えるでしょう。
以上で終わりにしたいと思います。このプレゼンテーションを聞いていただき、本当に嬉しく思います。ありがとうございました。
