知財求人 - 知財ポータルサイト「IP Force」
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2024077550
(43)【公開日】2024-06-07
(54)【発明の名称】シミュレーション装置及びシミュレーション方法
(51)【国際特許分類】
A61B 5/11 20060101AFI20240531BHJP
A61B 6/03 20060101ALI20240531BHJP
【ＦＩ】
A61B5/11 300
A61B6/03 360B
【審査請求】未請求
【請求項の数】9
【出願形態】ＯＬ
(21)【出願番号】P 2022189699
(22)【出願日】2022-11-28
(71)【出願人】
【識別番号】000006138
【氏名又は名称】株式会社明治
(71)【出願人】
【識別番号】512078812
【氏名又は名称】道脇 幸博
(74)【代理人】
【識別番号】110002675
【氏名又は名称】弁理士法人ドライト国際特許事務所
(72)【発明者】
【氏名】道脇 幸博
(72)【発明者】
【氏名】菊地 貴博
【テーマコード（参考）】
4C038
4C093
【Ｆターム（参考）】
4C038VA04
4C038VB08
4C038VB09
4C038VC15
4C093AA22
4C093AA26
4C093CA18
4C093DA04
4C093FF24
(57)【要約】
【課題】嚥下時又は咀嚼時における頭頸部器官の運動に関して、筋シナジーを含めたシミュレーションを行え得るシミュレーション装置及びシミュレーション方法を提供する。
【解決手段】嚥下シミュレーション装置１は、筋種ごとに筋粒子に筋活動率が設定された動的三次元頭頸部モデルを用いて、嚥下時における頭頸部器官の運動を粒子法に基づいてシミュレーションし、当該シミュレーション時における筋活動率の解析結果に基づいて嚥下時における頭頸部器官の運動に関して筋シナジーの解析を行え得る。
【選択図】図１
【解決手段】嚥下シミュレーション装置１は、筋種ごとに筋粒子に筋活動率が設定された動的三次元頭頸部モデルを用いて、嚥下時における頭頸部器官の運動を粒子法に基づいてシミュレーションし、当該シミュレーション時における筋活動率の解析結果に基づいて嚥下時における頭頸部器官の運動に関して筋シナジーの解析を行え得る。
【選択図】図１
【特許請求の範囲】
【請求項1】
複数の頭頸部器官が筋種ごとに三次元画像内で筋線維方向を特定した、粒子状又はメッシュ状の複数の筋形成体によりモデル化され、かつ、嚥下時又は咀嚼時における前記筋線維方向に基づく収縮応力が筋活動率により前記筋形成体に与えられる動的三次元頭頸部モデルを記憶する記憶部と、
前記筋活動率が変更可能な前記動的三次元頭頸部モデルを用い、嚥下時又は咀嚼時における前記頭頸部器官の運動を粒子法又は格子法に基づいて前記三次元画像でシミュレーションする運動解析部と、
前記頭頸部器官の運動の経時変化に応じて変化する各前記筋種の前記筋活動率から、前記筋種が協調運動する筋シナジーの個数と、嚥下時又は咀嚼時における前記筋シナジーごとの構成筋種の筋活動度合いとを示す空間パターンと、前記筋シナジーごとの活動度合いを経時変化で示した時間パターンと、を解析する筋シナジー解析処理部と、
を備える、シミュレーション装置。
前記筋活動率が変更可能な前記動的三次元頭頸部モデルを用い、嚥下時又は咀嚼時における前記頭頸部器官の運動を粒子法又は格子法に基づいて前記三次元画像でシミュレーションする運動解析部と、
前記頭頸部器官の運動の経時変化に応じて変化する各前記筋種の前記筋活動率から、前記筋種が協調運動する筋シナジーの個数と、嚥下時又は咀嚼時における前記筋シナジーごとの構成筋種の筋活動度合いとを示す空間パターンと、前記筋シナジーごとの活動度合いを経時変化で示した時間パターンと、を解析する筋シナジー解析処理部と、
を備える、シミュレーション装置。
【請求項2】
前記筋シナジー解析処理部は、
前記運動解析部による嚥下時又は咀嚼時における前記頭頸部器官の運動のシミュレーション結果から、嚥下時又は咀嚼時における前記筋種ごとの前記筋活動率の経時変化を表した筋活動行列Ｘを生成し、
設定された前記筋シナジーの個数を用いて前記筋活動行列Ｘを、前記筋シナジーごとに前記構成筋種ごとの前記筋活動度合いをそれぞれ表した空間パターン行列Ｗと、前記筋シナジーごとに活動度合いの経時変化を表した時間パターン行列Ｈとに変換し、前記空間パターン及び前記時間パターンの解析結果として、前記空間パターン行列Ｗと前記時間パターン行列Ｈとを算出する、請求項１に記載のシミュレーション装置。
前記運動解析部による嚥下時又は咀嚼時における前記頭頸部器官の運動のシミュレーション結果から、嚥下時又は咀嚼時における前記筋種ごとの前記筋活動率の経時変化を表した筋活動行列Ｘを生成し、
設定された前記筋シナジーの個数を用いて前記筋活動行列Ｘを、前記筋シナジーごとに前記構成筋種ごとの前記筋活動度合いをそれぞれ表した空間パターン行列Ｗと、前記筋シナジーごとに活動度合いの経時変化を表した時間パターン行列Ｈとに変換し、前記空間パターン及び前記時間パターンの解析結果として、前記空間パターン行列Ｗと前記時間パターン行列Ｈとを算出する、請求項１に記載のシミュレーション装置。
【請求項3】
前記筋シナジー解析処理部は、
嚥下時又は咀嚼時における前記頭頸部器官の特徴的な運動に基づいて、前記筋活動行列Ｘを時系列に沿って複数の時間領域に分割された複数の時分割筋活動行列Ｘｎを生成し、前記複数の時分割筋活動行列Ｘｎごとに、前記空間パターン行列Ｗ及び前記時間パターン行列Ｈを解析する、請求項２に記載のシミュレーション装置。
嚥下時又は咀嚼時における前記頭頸部器官の特徴的な運動に基づいて、前記筋活動行列Ｘを時系列に沿って複数の時間領域に分割された複数の時分割筋活動行列Ｘｎを生成し、前記複数の時分割筋活動行列Ｘｎごとに、前記空間パターン行列Ｗ及び前記時間パターン行列Ｈを解析する、請求項２に記載のシミュレーション装置。
【請求項4】
前記筋シナジー解析処理部は、
前記時間パターン行列Ｈを時系列に沿って複数の時間領域に分割して複数の時分割時間パターン行列Ｈｎを生成し、前記空間パターン行列Ｗが前記複数の時分割時間パターン行列Ｈｎに対応する構成を有する、請求項２に記載のシミュレーション装置。
前記時間パターン行列Ｈを時系列に沿って複数の時間領域に分割して複数の時分割時間パターン行列Ｈｎを生成し、前記空間パターン行列Ｗが前記複数の時分割時間パターン行列Ｈｎに対応する構成を有する、請求項２に記載のシミュレーション装置。
【請求項5】
前記筋シナジー解析処理部は、
前記時分割時間パターン行列Ｈｎの前記時間領域における前記筋種の前記筋活動率を補間する、請求項４に記載のシミュレーション装置。
前記時分割時間パターン行列Ｈｎの前記時間領域における前記筋種の前記筋活動率を補間する、請求項４に記載のシミュレーション装置。
【請求項6】
前記筋シナジー解析処理部は、
前記筋シナジーごとに得られた前記筋シナジーの活動度合いを集約した前記筋シナジーの解析結果を得、前記解析結果を反映した筋活動率を前記筋種ごとに生成し、
前記運動解析部は、
前記筋シナジー解析処理部で前記筋種ごとに得られた前記筋シナジーの解析結果を反映した筋活動率を前記筋種ごとに設定した前記動的三次元頭頸部モデルにより嚥下時又は咀嚼時における前記頭頸部器官の運動を前記粒子法又は前記格子法に基づいて前記三次元画像でシミュレーションする、請求項１に記載のシミュレーション装置。
前記筋シナジーごとに得られた前記筋シナジーの活動度合いを集約した前記筋シナジーの解析結果を得、前記解析結果を反映した筋活動率を前記筋種ごとに生成し、
前記運動解析部は、
前記筋シナジー解析処理部で前記筋種ごとに得られた前記筋シナジーの解析結果を反映した筋活動率を前記筋種ごとに設定した前記動的三次元頭頸部モデルにより嚥下時又は咀嚼時における前記頭頸部器官の運動を前記粒子法又は前記格子法に基づいて前記三次元画像でシミュレーションする、請求項１に記載のシミュレーション装置。
【請求項7】
前記筋シナジー解析処理部で前記筋種ごとに得られた前記筋シナジーの解析結果を反映した筋活動率を前記筋種ごとに設定した前記動的三次元頭頸部モデルにより、嚥下時又は咀嚼時における前記頭頸部器官の運動を前記粒子法又は前記格子法に基づいて前記三次元画像でシミュレーションした結果に基づいて、前記筋シナジーの個数を解析する、請求項６に記載のシミュレーション装置。
【請求項8】
前記運動解析部による嚥下時又は咀嚼時における前記頭頸部器官の運動のシミュレーション結果から、前記筋シナジー解析処理部による解析が必要であるか否かを判定するシミュレーション判定部を備える、請求項１に記載のシミュレーション装置。
【請求項9】
複数の頭頸部器官が筋種ごとに三次元画像内で筋線維方向を特定した、粒子状又はメッシュ状の複数の筋形成体によりモデル化され、かつ、嚥下時又は咀嚼時における前記筋線維方向に基づく収縮応力が筋活動率により前記筋形成体に与えられる動的三次元頭頸部モデルを記憶部に記憶する記憶ステップと、
運動解析部によって、前記筋活動率が変更可能な前記動的三次元頭頸部モデルを用い、嚥下時又は咀嚼時における前記頭頸部器官の運動を粒子法又は格子法に基づいて前記三次元画像でシミュレーションする運動解析ステップと、
前記頭頸部器官の運動の経時変化に応じて変化する各前記筋種の前記筋活動率から、前記筋種が協調運動する筋シナジーの個数と、嚥下時又は咀嚼時における前記筋シナジーごとの構成筋種の筋活動度合いを示す空間パターンと、前記筋シナジーごとの活動度合いを経時変化で示した時間パターンと、を筋シナジー解析処理部によって解析する筋シナジー解析処理ステップと、
を含む、シミュレーション方法。
運動解析部によって、前記筋活動率が変更可能な前記動的三次元頭頸部モデルを用い、嚥下時又は咀嚼時における前記頭頸部器官の運動を粒子法又は格子法に基づいて前記三次元画像でシミュレーションする運動解析ステップと、
前記頭頸部器官の運動の経時変化に応じて変化する各前記筋種の前記筋活動率から、前記筋種が協調運動する筋シナジーの個数と、嚥下時又は咀嚼時における前記筋シナジーごとの構成筋種の筋活動度合いを示す空間パターンと、前記筋シナジーごとの活動度合いを経時変化で示した時間パターンと、を筋シナジー解析処理部によって解析する筋シナジー解析処理ステップと、
を含む、シミュレーション方法。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、シミュレーション装置及びシミュレーション方法に関する。
【背景技術】
【0002】
嚥下時や咀嚼時の食品物性と頭頸部器官の運動との関係は複雑であり、頭頸部器官の運動を正確に把握することは非常に困難である。ここで、嚥下とは、口腔内に取り込まれた食品（飲料を含む）を、咽頭・食道を経て胃に送り込む運動である。嚥下時には、口腔、咽頭、喉頭、食道の各筋が、短時間のうちに決められた順序で活動し、複雑な運動を遂行している。また、咀嚼とは、口腔内に取り込まれた食品を、下顎骨やオトガイ舌筋等を動かしながら、かみ砕いてゆく運動である。咀嚼時においても、口腔、咽頭、喉頭の筋やオトガイ舌筋等の各筋が決められた順序で活動し、複雑な運動を遂行している。
【0003】
ところで、近年、このような筋活動解析として、被験者の大腿二頭筋や内転筋等の複数の部位にそれぞれ電極及び筋電計を貼付し、最大随意収縮（ＭＶＣ：Maximum Voluntary Contraction）時の筋電位を測定し、複数の筋種が協調運動によって制御されるとした筋シナジー仮説と、測定結果とに基づいて歩行時における筋シナジーの解析を行うことが考えられている（例えば、非特許文献１参照）。
【0004】
また、摂食嚥下の筋活動解析としては、舌骨上筋群用筋電センサ又は舌骨下筋群用筋電センサ（多チャンネルの電極が整列したアレイ状電極）を使用して、生体信号を検出し、筋シナジー仮説に基づいて変換された特徴量を用いることを特徴とした摂食嚥下機能評価方法が提案されている（例えば、特許文献１参照）。
【0005】
他方、摂食嚥下の筋活動解析としては、筋種ごとに筋活動率の値を設定する嚥下シミュレーション装置が提案されている（例えば、特許文献２参照）。
【先行技術文献】
【特許文献】
【0006】
【特許文献1】特開２０２１－１４２０８７号公報
【特許文献2】特開２０２１－１１２５５７号公報
【非特許文献】
【0007】
【非特許文献1】松永 直人，金岡 恒治，“早稲田大学審査学位論文 博士（スポーツ科学） シナジー解析を用いたスポーツ動作の筋活動解析”［令和4年11月18日検索］，インターネット＜URL：https://waseda.repo.nii.ac.jp/?action=repository_action_common_download&item_id=41273&item_no=1&attribute_id=20&file_no=1＞
【発明の概要】
【発明が解決しようとする課題】
【0008】
しかしながら、非特許文献１では、被検者の肌に筋電計等を直接貼付して筋電位等を測定することしかできないため、人体内部である口腔や咽頭、喉頭等の頭頸部器官については筋電位等を測定できず、嚥下時や咀嚼時における頭頸部器官の運動に関する筋シナジーの解析を行うことが困難であるという問題があった。
【0009】
特許文献１では、筋電センサを使用しても、舌、軟口蓋、咽頭の筋電位は取得できず、また、舌骨や喉頭が移動すること及びクロストーク問題などから、舌骨上筋群・下筋群の特定の筋電位は精密には特定できないと考えられるため、嚥下時や咀嚼時における頭頸部器官の運動に関する筋シナジーの解析を行うことが困難である。
【0010】
特許文献２では、筋種ごとに筋活動率の値を設定しており、筋シナジーは何ら考慮されていない。
【0011】
本発明は以上の点を考慮してなされたもので、嚥下時又は咀嚼時における頭頸部器官の運動に関して、筋シナジーを含めたシミュレーションを行え得るシミュレーション装置及びシミュレーション方法を提供することを目的とする。
【課題を解決するための手段】
【0012】
本発明に係るシミュレーション装置は、複数の頭頸部器官が筋種ごとに三次元画像内で筋線維方向を特定した、粒子状又はメッシュ状の複数の筋形成体によりモデル化され、かつ、嚥下時又は咀嚼時における前記筋線維方向に基づく収縮応力が筋活動率により前記筋形成体に与えられる動的三次元頭頸部モデルを記憶する記憶部と、前記筋活動率が変更可能な前記動的三次元頭頸部モデルを用い、嚥下時又は咀嚼時における前記頭頸部器官の運動を粒子法又は格子法に基づいて前記三次元画像でシミュレーションする運動解析部と、前記頭頸部器官の運動の経時変化に応じて変化する各前記筋種の前記筋活動率から、前記筋種が協調運動する筋シナジーの個数と、嚥下時又は咀嚼時における前記筋シナジーごとの構成筋種の筋活動度合いを示す空間パターンと、前記筋シナジーごとの活動度合いを経時変化で示した時間パターンと、を解析する筋シナジー解析処理部と、を備える。
【0013】
また、本発明に係るシミュレーション方法は、複数の頭頸部器官が筋種ごとに三次元画像内で筋線維方向を特定した、粒子状又はメッシュ状の複数の筋形成体によりモデル化され、かつ、嚥下時又は咀嚼時における前記筋線維方向に基づく収縮応力が筋活動率により前記筋形成体に与えられる動的三次元頭頸部モデルを記憶部に記憶する記憶ステップと、運動解析部によって、前記筋活動率が変更可能な前記動的三次元頭頸部モデルを用い、嚥下時又は咀嚼時における前記頭頸部器官の運動を粒子法又は格子法に基づいて前記三次元画像でシミュレーションする運動解析ステップと、前記頭頸部器官の運動の経時変化に応じて変化する各前記筋種の前記筋活動率から、前記筋種が協調運動する筋シナジーの個数と、嚥下時又は咀嚼時における前記筋シナジーごとの構成筋種の筋活動度合いを示す空間パターンと、前記筋シナジーごとの活動度合いを経時的変化で示した時間パターンと、を筋シナジー解析処理部によって解析する筋シナジー解析処理ステップと、を含む。
【発明の効果】
【0014】
本発明によれば、筋種ごとに筋形成体に設定された筋活動率を変更可能な動的三次元頭頸部モデルを用いて、嚥下時又は咀嚼時における頭頸部器官の運動に関して、筋シナジーを含めたシミュレーションを行え得る。
【図面の簡単な説明】
【0015】
【図1】嚥下シミュレーション装置の回路構成を示すブロック図である。
【図2】動的三次元頭頸部モデルの構成を示す概略図である。
【図3】動的三次元頭頸部モデルにおいて頭頸部器官を構成する粒子を概略的に示した概略図である。
【図5】５Ａは、動的三次元頭頸部モデルのオトガイ舌筋周辺の構成を示す概略図であり、５Ｂは、動的三次元頭頸部モデルの横舌筋の構成を示す概略図である。
【図6】動的三次元頭頸部モデルの顎舌骨筋、舌骨及び下顎骨の構成を示す概略図である。
【図7】動的三次元頭頸部モデルの口蓋咽頭筋周辺の構成を示す概略図である。
【図8】８Ａは、動的三次元頭頸部モデルの舌骨、甲状軟骨及び輪状軟骨、下咽頭収縮筋甲状咽頭部周辺の構成を示す概略図であり、８Ｂは、動的三次元頭頸部モデルの中咽頭収縮筋大角咽頭部の構成を示す概略図であり、８Ｃは、動的三次元頭頸部モデルの中咽頭収縮筋小角咽頭部の構成を示す概略図である。
【図9】嚥下時における動的三次元頭頸部モデルの状態変化を示した概略図である。
【図10】筋粒子について説明するための概略図である。
【図11】筋粒子の筋線維方向を説明するための概略図である。
【図12A】動的三次元頭頸部モデルの各頭頸部器官の嚥下時における筋活動率の経時変化を示すグラフ（１）である。
【図12B】動的三次元頭頸部モデルの各頭頸部器官の嚥下時における筋活動率の経時変化を示すグラフ（２）である。
【図13】非負値行列因子分解を説明するための概略図である。
【図14】第１の筋シナジー、第２の筋シナジー及び第３の筋シナジーの構成筋種と、算出された筋種ごとの重みとを示す概略図である。
【図15】第４の筋シナジー、第５の筋シナジー及び第６の筋シナジーの構成筋種と、算出された筋種ごとの重みとを示す概略図である。
【図16】第７の筋シナジー、第８の筋シナジー及び第９の筋シナジーの構成筋種と、算出された筋種ごとの重みとを示す概略図である。
【図17】筋シナジーの活動度合いの経時変化である時間パターンを示すグラフである。
【図18】１８Ａは、複数の時間領域に分割した筋活動行列Ｘ及び時間パターン行列Ｈを説明するための概略図であり、１８Ｂは、筋シナジーの空間パターン行列Ｗ及び時間パターン行列Ｈの構成を説明するための概略図である。
【図19】動的三次元頭頸部モデルでの舌骨の運動とその変位との関係を示すとともに、筋活動行列Ｘを分割した分割筋活動行列Ｘ１～６の時間領域を示したグラフである。
【図20】２０Ａは、動的三次元頭頸部モデルの嚥下時におけるオトガイ舌筋等の舌を構成する筋について筋活動率の経時変化を示すグラフであり、２０Ｂは、動的三次元頭頸部モデルの嚥下時における口蓋帆挙筋等の軟口蓋を構成する筋について筋活動率の経時変化を示すグラフであり、２０Ｃは、動的三次元頭頸部モデルの嚥下時における顎二腹筋前腹等の舌骨上筋群を構成する筋について筋活動率の経時変化を示すグラフであり、２０Ｄは、動的三次元頭頸部モデルの嚥下時における口蓋咽頭筋縦走部等の咽頭を構成する筋について筋活動率の経時変化を示すグラフである。
【図21】２１Ａは、筋活動行列Ｘを時間分割した時分割筋活動行列Ｘ１と空間パターン行列Ｗ１と時分割時間パターン行列Ｈ１との関係を説明するための概略図であり、２１Ｂは、筋活動行列Ｘを時間分割した時分割筋活動行列Ｘ２と空間パターン行列Ｗ２と時分割時間パターン行列Ｈ２との関係を説明するための概略図であり、２１Ｃは、筋活動行列Ｘを時間分割した時分割筋活動行列Ｘ４と空間パターン行列Ｗ４～７と時分割時間パターン行列Ｈ４との関係を説明するための概略図である。
【図22】筋活動率を補間する補間工程を説明するための表である。
【発明を実施するための形態】
【0016】
以下図面について本発明の一実施の形態を詳述する。以下の説明において、同一の構成要素には同一の符号を付し、重複する説明は省略する。
【0017】
（１）＜本発明の概略について＞
図１は本実施形態に係る嚥下シミュレーション装置１の回路構成を示したブロック図である。嚥下シミュレーション装置１は、パーソナルコンピュータ（ＰＣとも称する）２と、入力部８１と、表示部４と、記憶部８３とを備えている。入力部８１は、マウス、キーボード等の入力機器であり、開発者からの操作命令をパーソナルコンピュータ２に対して入力し、パーソナルコンピュータ２において操作命令に応じた各種演算処理を実行させる。記憶部８３は、パーソナルコンピュータ２にて形成した粒子表示の動的三次元頭頸部モデル（図２等において後述する）や、各種設定条件、解析結果、後述する筋活動率等を記憶する。
図１は本実施形態に係る嚥下シミュレーション装置１の回路構成を示したブロック図である。嚥下シミュレーション装置１は、パーソナルコンピュータ（ＰＣとも称する）２と、入力部８１と、表示部４と、記憶部８３とを備えている。入力部８１は、マウス、キーボード等の入力機器であり、開発者からの操作命令をパーソナルコンピュータ２に対して入力し、パーソナルコンピュータ２において操作命令に応じた各種演算処理を実行させる。記憶部８３は、パーソナルコンピュータ２にて形成した粒子表示の動的三次元頭頸部モデル（図２等において後述する）や、各種設定条件、解析結果、後述する筋活動率等を記憶する。
【0018】
パーソナルコンピュータ２は、例えば、頭頸部器官からなる動的三次元頭頸部モデル（図２等において後述する）を三次元画像により形成し、経口摂取品を三次元画像内で擬似経口摂取品としてモデル化する。パーソナルコンピュータ２は、動的三次元頭頸部モデルにおける各頭頸部器官の運動と、擬似経口摂取品の嚥下時の挙動とを、粒子法を用いて三次元画像内で解析する。
【0019】
なお、本実施形態では、擬似経口摂取品の嚥下時の頭頸部器官の挙動を動的三次元頭頸部モデルで再現して嚥下時における筋シナジーの解析を行う場合について以下説明するが、本発明はこれに限らず、擬似経口摂取品の咀嚼時の頭頸部器官の挙動を動的三次元頭頸部モデルで再現して咀嚼時における筋シナジーの解析を行うようにしてもよい。
【0020】
また、パーソナルコンピュータ２は、擬似経口摂取品を三次元画像内に設けずに、動的三次元頭頸部モデルにおける嚥下時の各頭頸部器官の運動のみを粒子法を用いて三次元画像内で解析することもできる。
【0021】
このような動的三次元頭頸部モデルは、例えば、誤嚥をし易い嚥下障害者の頭頸部、又は、誤嚥をし難い健常者の頭頸部等を模倣して形成する。例えば、誤嚥をし易い嚥下障害者の頭頸部をモデル化した動的三次元頭頸部モデルでは、嚥下障害者の嚥下時における各頭頸部器官の運動や、擬似経口摂取品の嚥下時の挙動を解析し、嚥下障害者に関して嚥下時の筋シナジーを含めた嚥下シミュレーションを行うことができる。一方、誤嚥をし難い健常者の頭頸部をモデル化した動的三次元頭頸部モデルでは、健常者の嚥下時における各頭頸部器官の運動や、擬似経口摂取品の嚥下時の挙動を解析し、健常者に関して嚥下時の筋シナジーを含めた嚥下シミュレーションを行うことができる。
【0022】
パーソナルコンピュータ２で得られた解析結果は、表示部４に出力され、表示部４の表示画面に表示される。表示部４は、例えばディスプレイ等であり、パーソナルコンピュータ２から出力された動的三次元頭頸部モデルの三次元画像や、擬似経口摂取品、解析結果、筋種ごとの筋活動率、複数の筋種が協調運動によって制御されるとした筋シナジー仮説に基づいて特定した複数の筋シナジーの空間パターンと時間パターン等を表示画面に表示する。空間パターンと時間パターンとについては後述する。
【0023】
これにより、表示部４は、動的三次元頭頸部モデルにおける各頭頸部器官の運動や、擬似経口摂取品の嚥下時の挙動、筋種ごとの筋活動率、筋シナジーに関する各種解析結果等を、開発者に対し視認させることができる。なお、本実施形態で三次元画像とは、動的三次元頭頸部モデルを仮想三次元空間に表現した動画像、及び、動画像を構成するフレーム画像を示す。
【0024】
嚥下シミュレーション装置１では、擬似経口摂取品の食塊量や粘度、比重等の物性値を変えて、動的三次元頭頸部モデルによる嚥下シミュレーションを行うことができ、動的三次元頭頸部モデルにより誤嚥の有無等を確認することや、筋シナジーを解析することができる。また、嚥下シミュレーション装置１では、三次元画像内に擬似経口摂取品を設定しない場合でも、動的三次元頭頸部モデルによる嚥下シミュレーションを行うことができ、シミュレーション結果に基づいて嚥下時における動的三次元頭頸部モデルの各頭頸部器官の運動に関して筋シナジーを解析することができる。
【0025】
本実施形態では、液面の変形や飛沫等の表現が容易な解析方法として、解析対象の液体や固体を粒子として扱う粒子法を用い、この粒子法によって、動的三次元頭頸部モデルにおける頭頸部器官の運動（動作）や、経口摂取品の挙動を、三次元画像内に表して嚥下シミュレーションを行なう。粒子法としては、特にＭＰＳ（Moving Particle Semi-implicit）法（Koshizuka et al，Comput．Fluid Dynamics J，4，29-46，1995）を適用することが望ましい。嚥下シミュレーションによって嚥下時における擬似経口摂取品の挙動を解析する際の粒子法としては、ＭＰＳ（Moving Particle Semi-implicit）法、ハミルトニアン粒子法（Hamiltonian MPS法：HMPS法）、又は離散化要素法（Discrete Element Method）を適用することが望ましい。
【0026】
また、嚥下シミュレーションによって動的三次元頭頸部モデルにおける各粒子の運動を解析する際の粒子法としては、ハミルトニアン粒子法（Hamiltonian MPS法）を適用することが望ましい。本実施形態では、嚥下シミュレーションによって動的三次元頭頸部モデルにおける各粒子の運動を解析する粒子法として、ハミルトニアン粒子法（Hamiltonian MPS法）を適用した場合について以下説明する。
【0027】
本実施形態の粒子法では、擬似経口摂取品を粒子に置き換えるだけでなく、動的三次元頭頸部モデルにおける頭頸部器官についても粒子に置き換え、粒子ごとに物理量を計算する。その結果、動的三次元頭頸部モデルにおける頭頸部器官や、擬似経口摂取品の嚥下時における微妙な変化の解析が可能となる。なお、本実施形態では、このような動的三次元頭頸部モデルを用いて行われる、粒子法による演算処理を、嚥下シミュレーションとも称する。
【0028】
ここで、本実施形態に係る嚥下シミュレーション装置１では、動的三次元頭頸部モデルにおける頭頸部器官について粒子に置き換えるだけでなく、さらに、口腔、咽頭部、喉頭部等の頭頸部器官を形成する粒子に対して、医学的知見に基づいて頭頸部器官の筋種ごとに三次元画像内で筋線維方向が特定され、かつ筋線維方向に基づく収縮応力が与えられている。
【0029】
本実施形態では、このように頭頸部器官の筋種ごとに三次元画像内で筋線維方向が特定され、かつ筋線維方向に基づく収縮応力が与えられる粒子を筋粒子（粒子状の筋形成体を筋粒子）と称する。そして、嚥下シミュレーション装置１では、筋種ごとに設定した筋粒子の嚥下時における挙動（筋線維方向に与えられる収縮応力の時間的変化）を、時系列に変化する筋活動率の値を調整することで再現している。
【0030】
ここで、筋活動率は、嚥下時に筋線維方向に基づく収縮応力が与えられる筋粒子に関し、筋種ごとに設定される、嚥下時における筋粒子の収縮応力の時間的変化を示すものである。嚥下シミュレーション装置１では、嚥下時に時系列に変化する筋活動率の値を筋種ごとにそれぞれ設定することで、筋種ごとに嚥下時における筋粒子の挙動を変えることができる。従って、筋種ごとに筋活動率の値を最適な値に調整することで、嚥下時における各頭頸部器官の理想的な動きを実現することができる。
【0031】
嚥下シミュレーション装置１では、所定の値の筋活動率を筋種ごとに設定した動的三次元頭頸部モデルを用いて嚥下シミュレーション（粒子法による演算処理）を行い、嚥下時の筋粒子の挙動に関する解析結果（筋粒子の位置と速度）を得る。ここで、本実施形態に係る嚥下シミュレーション装置１にて使用する動的三次元頭頸部モデルは、例えば、特開２０２１－０２９９８０号公報に開示された手法に従って動的三次元頭頸部モデルを作製し、当該動的三次元頭頸部モデルを形成する粒子に対して筋種ごとにそれぞれ筋活動率を設定した動的三次元頭頸部モデルを用いることができる。また、本実施形態に係る嚥下シミュレーション装置１にて使用する動的三次元頭頸部モデルは、特開２０２１－１１２５５７号公報で開示された手法に従って筋活動率が設定された動的三次元頭頸部モデルを用いることができる。
【0032】
本実施形態において、動的三次元頭頸部モデルを形成する粒子に筋種ごとにそれぞれ設定される筋活動率は、嚥下シミュレーション時に動的三次元頭頸部モデルにより嚥下時の動作を再現するように開発者が手動で筋種ごとに筋活動率を設定してもよく、また、上記の特開２０２１－１１２５５７号公報で開示されているように評価関数の解析結果を基に、筋種ごとに最適な値に筋活動率の値を設定するようにしてもよい。
【0033】
嚥下シミュレーション装置１では、筋種ごとに筋活動率の値を適宜調整することにより、嚥下時における各頭頸部器官の理想的な動きを実現し得る。そして、嚥下時における頭頸部器官の運動や、擬似経口摂取品の挙動を正確に再現した際の筋種ごとの筋活動率から、最適な筋シナジーの個数や、各筋シナジーの空間パターンと時間パターンとを特定することができる。
【0034】
嚥下シミュレーション装置１では、得られた筋シナジーの空間パターンと時間パターンを基に各筋種の筋活動率を修正し、再度、動的三次元頭頸部モデルを動作させる。これにより、動的三次元頭頸部モデルの動作が所望する嚥下動作となっている場合には、筋シナジーによる個数が最適であると推認し得、嚥下時に協調運動によって制御される筋種の組み合わせ数を筋シナジーの個数から特定することができる。
【0035】
（２）＜筋シナジー解析装置におけるパーソナルコンピュータの回路構成＞
次に、嚥下シミュレーション装置１のパーソナルコンピュータ２について以下説明する。図１に示すように、パーソナルコンピュータ２は、頭頸部モデリング部１０、器官運動設定部３０、経口摂取品物性設定部４０、運動解析部５０、物性特定部７０、制御部９０及び筋活動率変更部８０、筋活動率解析部９２、筋シナジー解析処理部９３及びシミュレーション判定部９４を備えている。頭頸部モデリング部１０は、例えば、図２に示すように、頭頸部器官を忠実に再現した動的三次元頭頸部モデル１０ｃを三次元画像により形成する。
次に、嚥下シミュレーション装置１のパーソナルコンピュータ２について以下説明する。図１に示すように、パーソナルコンピュータ２は、頭頸部モデリング部１０、器官運動設定部３０、経口摂取品物性設定部４０、運動解析部５０、物性特定部７０、制御部９０及び筋活動率変更部８０、筋活動率解析部９２、筋シナジー解析処理部９３及びシミュレーション判定部９４を備えている。頭頸部モデリング部１０は、例えば、図２に示すように、頭頸部器官を忠実に再現した動的三次元頭頸部モデル１０ｃを三次元画像により形成する。
【0036】
ここでは、図２に示す動的三次元頭頸部モデル１０ｃは、例えば、図３に示すように、粒子表示の動的三次元頭頸部モデル１０ｃからマーチングキューブ法等を用いて作製された構成を有する。このような動的三次元頭頸部モデル１０ｃは、当該動的三次元頭頸部モデル１０ｃに設定した各粒子の設定情報はそのままに、粒子１つ１つを表示せずに、単に頭頸部器官の表面を表示して簡易的に表している。
【0037】
また、図２に示すように、三次元画像内に形成する擬似経口摂取品１００についても、表示部４に表示させる際には、擬似経口摂取品１００を形成している粒子１つ１つは表示させずに、マーチングキューブ法等により生成した、擬似経口摂取品１００の表面形状のみを表示させる。
【0038】
なお、本実施形態においては、粒子１つ１つを表示せずに単に表面を表示して簡易的に表した動的三次元頭頸部モデル１０ｃや擬似経口摂取品１００を適用した例を用いて説明するが、本発明はこれに限らず、例えば、図３に示すように、粒子１つ１つを表示した動的三次元頭頸部モデル１０ｃや擬似経口摂取品１００を適用してもよい。
【0039】
ここで、図１に示す器官運動設定部３０は、動的三次元頭頸部モデル１０ｃにおける各頭頸部器官の運動を設定する。本実施形態における器官運動設定部３０は、強制運動設定部３１と、筋収縮運動設定部３２とを備える。強制運動設定部３１は、動的三次元頭頸部モデル１０ｃによる擬似経口摂取品１００の嚥下時に、頭頸部器官で強制的に移動する粒子を強制移動粒子として設定し、これら複数の強制移動粒子の運動を設定する。筋収縮運動設定部３２は、医学的知見に基づき頭頸部器官の筋種ごとに三次元画像内で筋線維方向が特定され、かつ筋線維方向に基づく収縮応力が与えられる粒子を筋粒子として設定し、擬似経口摂取品１００の嚥下時における当該収縮応力に基づいて筋粒子の運動を設定する。
【0040】
本実施形態に係る動的三次元頭頸部モデル１０ｃを形成する粒子は、強制移動粒子と、筋粒子と、これら強制移動粒子及び筋粒子以外の粒子と、の３種類のいずれかに定義される。動的三次元頭頸部モデル１０ｃは、器官運動設定部３０による設定状態を基に、各頭頸部器官が動いて嚥下シミュレーションを実行することができる。なお、強制移動粒子と、筋粒子と、これら強制移動粒子及び筋粒子以外の粒子とについて、特に区別する必要がない場合には、以下、強制移動粒子と、筋粒子と、これら強制移動粒子及び筋粒子以外の粒子と、の３種の粒子をまとめて、単に粒子と称する。
【0041】
経口摂取品物性設定部４０は、解析対象としての飲食品、医薬品又は医薬部外品等の経口摂取品の物性値を設定し、経口摂取品をモデル化した擬似経口摂取品１００を三次元画像内に形成する。なお、経口摂取品物性設定部４０は、解析対象として異なる物性の液体、半固体又は固体の複数の擬似経口摂取品１００を設定することができる。なお、半固体としては例えばゼリー等を含み、固体としては例えば錠剤等も含む。
【0042】
本実施形態の場合、経口摂取品物性設定部４０は、経口摂取品の物性値として、例えば、経口摂取品となる食塊の密度［ｇ／ｍＬ］と、動的三次元頭頸部モデル１０ｃに嚥下させる食塊量［ｍＬ］と、表面張力［Ｎ／ｍ］と、各頭頸部器官における接触角と、各頭頸部器官におけるスリップ係数と、を設定する。なお、ここでスリップ係数とは、生体表面と食塊（経口摂取品）の表面の濡れ性、撥水性を制御するパラメータであり、接触面における見かけの粘度として考えることができる。スリップ係数が大きい場合は界面での摩擦が大きくなり、結果的に食塊の動きにブレーキをかける効果がある。スリップ係数が小さい場合は界面での摩擦が小さくなり、０の場合は鏡面のような状態となる。スリップ係数１は流体の粘度と同等程度の摩擦効果を界面に与えることを意味する。スリップ係数は、想定する経口摂取品が有する濡れ性や撥水性等を解析して決定する。
【0043】
経口摂取品物性設定部４０は、各頭頸部器官における接触角として、動的三次元頭頸部モデル１０ｃにおける咽頭、喉頭、舌、軟口蓋等での接触角をそれぞれ設定する。また、経口摂取品物性設定部４０は、各頭頸部器官におけるスリップ係数として、動的三次元頭頸部モデル１０ｃの咽頭、喉頭、舌、軟口蓋等でのスリップ係数をそれぞれ設定する。
【0044】
なお、本実施形態においては、擬似経口摂取品１００の物性値として上述した物性値のみだけでなく、例えば、擬似経口摂取品１００が液体のときは、液量・粘度・表面張力・比重・熱伝導率・比熱等の物性値を設定するようにしてもよい。また、擬似経口摂取品１００が固体のときには、形状・寸法・弾性係数・引っ張り強さ・降伏点・降伏応力・粘度のずり速度依存性・動的粘弾性・静的粘弾性・圧縮応力・破断応力・破断ひずみ・硬度・付着性・凝集性・熱伝導率・比熱等の物性値を設定するようにしてもよい。さらに、擬似経口摂取品１００が半固体（可塑性があるが、流動性はない）であるときには、量・粘度・比重・降伏点・降伏点応力・粘度のずり速度依存性・動的粘弾性・静的粘弾性・圧縮応力・付着性・凝集性等の物性値を設定するようにしてもよい。
【0045】
運動解析部５０では、粒子法を用いた嚥下シミュレーションを行い、動的三次元頭頸部モデル１０ｃの嚥下時における頭頸部器官の運動と、頭頸部器官の運動に伴う擬似経口摂取品１００の嚥下時の挙動と、を解析する。
【0046】
ここで、図３は、図２に示した動的三次元頭頸部モデル１０ｃを構成する頭頸部器官についてマーチングキューブ法による処理をしないで粒子１つ１つを表示させた動的三次元頭頸部モデル１０ｃの一例を示した概略図であり、図４は、図３に示した動的三次元頭頸部モデル１０ｃの正中面における断面構成を示した断面図である。図２、図３及び図４に示す動的三次元頭頸部モデル１０ｃにおいて、粒子法による解析によって、舌１２の進行波的波動運動、喉頭蓋１５ａの回転運動、喉頭１５の往復運動、咽頭部１４の筋収縮運動等の動きがそれぞれ再現され、口腔から頭頸部内部に投入された擬似経口摂取品１００（図２及び図３では図示せず）を動かす。擬似経口摂取品１００の動きも粒子法により解析される。擬似経口摂取品１００は固体・半固体・液体のいずれでも粒子として取り扱われる。
【0047】
運動解析部５０は、例えば、経口摂取品物性設定部４０による擬似経口摂取品１００の物性値の変更や、筋活動率変更部８０による筋活動率の値の変更等が行われた後、粒子法を用いた嚥下シミュレーションを行う。運動解析部５０は、擬似経口摂取品１００の物性値が変更されると、粒子法を用いた嚥下シミュレーションによって、擬似経口摂取品１００が嚥下される際の経路が変化した解析結果を得ることができる。
【0048】
また、運動解析部５０は、例えば、筋活動率変更部８０により筋種ごとに筋活動率の値が変更されると、粒子法を用いた嚥下シミュレーションによって、舌１２の進行波的波動運動や、軟口蓋１３ｂの挙上運動、喉頭蓋１５ａの反転運動、喉頭１５の挙上運動、声帯１５ｃの内転運動、披裂部１５ｂの前方運動、咽頭部１４の収縮と挙上運動等の挙動がそれぞれ変化した解析結果を得ることができる。
【0049】
物性特定部７０は、運動解析部５０の解析結果を基に、誤嚥を回避できる擬似経口摂取品の食塊量、粘度及びせん断速度を推測する。このうち擬似経口摂取品の食塊量及び粘度は、経口摂取品物性設定部４０により設定される物性値である。
【0050】
制御部９０は、パーソナルコンピュータ２の各部を制御して、嚥下シミュレーション装置１の諸機能を実行させる。制御部９０は内蔵メモリに嚥下シミュレーター（解析用ソフトウェア）を保有する。
【0051】
筋活動率変更部８０は、入力部８１からの操作命令に従って、動的三次元頭頸部モデル１０ｃにおいて、筋種ごとに設定される筋活動率の値を変更する。これにより、運動解析部５０は、筋種ごとに筋活動率の値が変更された動的三次元頭頸部モデル１０ｃにより、粒子法を用いた嚥下シミュレーションを行える。
【0052】
筋活動率解析部９２は、例えば、運動解析部５０による動的三次元頭頸部モデル１０ｃの嚥下シミュレーションにより得られた解析結果と、そのときに設定されていた筋活動率の値とを用いて、当該嚥下シミュレーション時に筋種ごとに設定された筋活動率の値を評価する評価関数を算出し、嚥下シミュレーションごとに算出した複数の評価関数について、評価関数の値がほぼ変わらない定常状態になったか否かの判定等を行う。筋活動率解析部９２は、得られた評価関数及び筋活動率を利用して最適な筋活動率を探索する。なお、筋活動率解析部９２による筋活動率の解析手法は、既に、特開２０２１－１１２５５７号公報で公知であるため、ここでは詳細な説明を省略する。
【0053】
筋シナジー解析処理部９３は、運動解析部５０による嚥下時における動的三次元頭頸部モデル１０ｃの各頭頸部器官の運動のシミュレーション結果を基に、当該頭頸部器官の運動の経時変化に応じて変化する各筋種の筋活動率を取得する。筋シナジー解析処理部９３は、筋種ごとに取得した筋活動率から、複数の筋種が協調運動によって制御されるとした筋シナジー仮説に基づいて、筋種が協調運動する筋シナジーの個数を特定可能な構成を有する。この場合、筋シナジー解析処理部９３は、運動解析部５０による嚥下時における動的三次元頭頸部モデル１０ｃの頭頸部器官の運動のシミュレーション結果から、嚥下時における筋種ごとの筋活動率の経時変化を表した筋活動行列Ｘを生成する。
【0054】
開発者は、運動解析部５０によって得られた動的三次元頭頸部モデル１０ｃの嚥下シミュレーションや、筋シナジー解析処理部９３で得られた筋活動行列Ｘ等を確認し、複数の筋種が協調運動によって制御されるとした筋シナジー仮説に基づいて、協調運動する筋種の組み合わせ数を推測する。筋シナジー解析処理部９３には、開発者が推測した筋シナジーの個数が設定される。
【0055】
なお、本実施形態においては、運動解析部５０によって得られた動的三次元頭頸部モデル１０ｃの嚥下シミュレーションや、筋シナジー解析処理部９３で得られた筋活動行列Ｘ等を開発者が確認して筋シナジーの個数を推測するようにした場合について説明するが、本発明はこれに限らない。例えば、運動解析部５０によって得られた動的三次元頭頸部モデル１０ｃの嚥下シミュレーションや、筋シナジー解析処理部９３で得られた筋活動行列Ｘに対する筋シナジーの個数のパターンを、機械学習モデル等により学習させておき、これら嚥下シミュレーションや筋活動行列Ｘから自動的に筋シナジーの個数を定めるようにしてもよい。
【0056】
筋シナジー解析処理部９３は、設定された筋シナジーの個数に基づいて筋活動行列Ｘに対して非負値行列因子分解(ＮＭＦ：Nonnegative Matrix Factorization)を実行することにより、筋活動行列Ｘを空間パターン行列Ｗと時間パターン行列Ｈとに変換する。ここで、空間パターン行列Ｗは、筋シナジーごとの構成筋種の筋活動度合いをそれぞれ表した行列状の数値データであり、時間パターン行列Ｈは、筋シナジーごとの活動度合いの経時変化のデータを表した行列状の数値データである。
【0057】
筋シナジー解析処理部９３は、筋活動行列Ｘを、空間パターン行列Ｗと時間パターン行列Ｈの積で近似して表すように演算処理を行い、得られた空間パターン行列Ｗと時間パターン行列Ｈとを表示部４に表示させ、開発者へ提示する。これにより開発者は、提示された空間パターン行列Ｗと時間パターン行列Ｈと嚥下シミュレーションによる動的三次元頭頸部モデル１０ｃの動作状況とに基づいて、筋シナジーの個数や、各筋シナジーの空間パターン行列Ｗ、時間パターン行列Ｈを基に、嚥下時における頭頸部器官の運動を解析することができる。
【0058】
また、シミュレーション判定部９４は、運動解析部５０による嚥下時における動的三次元頭頸部モデル１０ｃの頭頸部器官の運動を示す嚥下シミュレーション結果から、筋シナジー解析処理部９３による解析が必要であるか否かを判定する。具体的には、運動解析部５０によって得られた動的三次元頭頸部モデル１０ｃの嚥下シミュレーションが最適な嚥下動作であるか否かを機械学習モデル等により学習させておき、機械学習モデル等に基づいて嚥下シミュレーション結果から自動的に筋シナジーの個数の変更が必要か否か等を判定する。例えば、運動解析部５０によって得られた動的三次元頭頸部モデル１０ｃの嚥下シミュレーションが最適な嚥下動作でないときには、筋シナジー解析処理部９３において設定する筋シナジーの個数の修正等を促し、修正した個数によって筋活動行列Ｘに対して非負値行列因子分解を実行させるように促すことができる。
【0059】
（３）＜動的三次元頭頸部モデルの構成＞
ここでは、本実施形態の特徴的構成である筋シナジー解析処理部９３を説明する前に、始めに、筋シナジーの解析に用いる動的三次元頭頸部モデル１０ｃについて以下説明する。
ここでは、本実施形態の特徴的構成である筋シナジー解析処理部９３を説明する前に、始めに、筋シナジーの解析に用いる動的三次元頭頸部モデル１０ｃについて以下説明する。
【0060】
上述したように、表示部４には、例えば、図２に示すような動的三次元頭頸部モデル１０ｃ及び擬似経口摂取品１００が表示され、運動解析部５０による解析結果として、動的三次元頭頸部モデル１０ｃにおいて擬似経口摂取品１００を嚥下したときの頭頸部器官の運動と、擬似経口摂取品１００の挙動と、が動画像により表示される。
【0061】
なお、本実施形態では、図２に示すように、粒子表示がされていない動的三次元頭頸部モデル１０ｃ及び擬似経口摂取品１００を用いて、頭頸部器官の運動と、擬似経口摂取品１００の挙動と、を動画像により提示する場合について説明するが、本発明はこれに限らず、図２及び図３に示すように、粒子１つ１つが表示された動的三次元頭頸部モデル１０ｃ及び擬似経口摂取品１００を用いて、頭頸部器官の運動と、擬似経口摂取品１００の挙動と、を動画像により提示するようにしてもよい。
【0062】
図２、図３及び図４において、Ｘ軸は、動的三次元頭頸部モデル１０ｃの正中面と直交する身体左右方向（正中面の面法線方向）を示し、Ｙ軸は、Ｘ軸と直交し、かつ動的三次元頭頸部モデル１０ｃの正中面と平行な身体前後方向を示し、Ｚ軸は、Ｘ軸及びＹ軸と直交する身体上下方向を示す。
【0063】
図２、図３及び図４に示すように、動的三次元頭頸部モデル１０ｃは、頭頸部器官として、甲状軟骨１１ｂ（図４）、オトガイ舌骨筋及びオトガイ舌筋等を含む舌１２、喉頭１５、声帯１５ｃ（図４）、披裂部１５ｂ（図４）、喉頭蓋１５ａ、気管１６（図４）、咽頭部１４（咽頭の管壁１８ａ、咽頭の粘膜１８ｂを含む）、口蓋１３（硬口蓋１３ａ及び軟口蓋１３ｂを含む）、食道１７（食道入口部１７ａ、食道の管壁１７ｂを含む）、下顎骨２１ａ等を有している。
【0064】
また、動的三次元頭頸部モデル１０ｃは、図５の５Ａに示すように、オトガイ舌筋がオトガイ舌筋（１）～（６）の６つの扇状領域に分割され、さらに、図５の５Ｂに示すように、横舌筋が身体前後方向で横舌筋（１）及び横舌筋（２）の２領域に分割されており、領域ごとにこれら筋種を構成する粒子に筋活動率がそれぞれ個別に設定可能に構成されている。さらに、動的三次元頭頸部モデル１０ｃには、その他、図５Ａに示すように、甲状軟骨が設けられているとともに、図６に示すように、舌骨や下顎骨等も設けられている。
【0065】
また、動的三次元頭頸部モデル１０ｃは、図６に示すように、顎舌骨筋が顎舌骨筋（１）及び顎舌骨筋（２）の２領域に分割され、当該領域ごとに異なる筋活動率が設定可能な構成を有している。さらに、動的三次元頭頸部モデル１０ｃには、図７に示すように、口蓋咽頭筋縦走部が口蓋咽頭筋縦走部（１）及び口蓋咽頭筋縦走部（２）の２領域に分割され、当該領域ごとに異なる筋活動率が設定可能な構成を有し、また、口蓋帆張筋、翼突鈎及び口蓋帆挙筋も設けられている。
【0066】
動的三次元頭頸部モデル１０ｃは、図８の８Ａ、８Ｂ及び８Ｃに示すように、舌骨及び甲状軟骨に加えて、輪状軟骨、下咽頭収縮筋甲状咽頭部、中咽頭収縮筋大角咽頭部及び中咽頭収縮筋小角咽頭部が設けられている。下咽頭収縮筋甲状咽頭部は、図８の８Ａに示すように、下咽頭収縮筋甲状咽頭部（１）及び下咽頭収縮筋甲状咽頭部（２）の２領域に分割され、中咽頭収縮筋大角咽頭部は、図８の８Ｂに示すように、中咽頭収縮筋大角咽頭部（１）及び中咽頭収縮筋大角咽頭部（２）の２領域に分割され、中咽頭収縮筋小角咽頭部は、図８の８Ｃに示すように、中咽頭収縮筋小角咽頭部（１）及び中咽頭収縮筋小角咽頭部（２）の２領域に分割されており、これら筋種ごとにそれぞれ異なる筋活動率が設定可能な構成を有する。
【0067】
本実施形態では、主に、舌１２、喉頭１５、声帯１５ｃ、披裂部１５ｂ、喉頭蓋１５ａ、気管１６、咽頭部１４、口蓋１３、食道１７、オトガイ舌筋、舌骨舌筋、下縦舌筋、茎突舌筋、上縦舌筋、横舌筋、垂直舌筋、口蓋帆挙筋、口蓋舌筋、口蓋帆張筋、口蓋垂筋、顎二腹筋前腹、オトガイ舌骨筋、顎舌骨筋、顎二腹筋後腹、茎突舌骨筋、胸骨舌骨筋、胸骨甲状筋、甲状舌骨筋、口蓋咽頭筋縦走部、口蓋咽頭筋横走部、茎突咽頭筋、上咽頭収縮筋翼突・頬咽頭部、上咽頭収縮筋舌咽頭部、中咽頭収縮筋小角咽頭部、中咽頭収縮筋大角咽頭部、下咽頭収縮筋甲状咽頭部、下咽頭収縮筋輪状咽頭部、披裂喉頭蓋筋等をまとめて頭頸部器官と称するが、これら１つ１つについても単に頭頸部器官とも称する。
【0068】
本実施形態では、擬似経口摂取品１００を粒子により表現するとともに、動的三次元頭頸部モデル１０ｃにおける頭頸部器官を粒子により表現する。ただし、上述したように、動的三次元頭頸部モデル１０ｃによる擬似経口摂取品１００の嚥下シミュレーションの解析結果を、動画像により開発者等に視認させる際には、動的三次元頭頸部モデル１０ｃの頭頸部器官と、擬似経口摂取品１００との表面形状をマーチングキューブ法等により表示させることが望ましい。これにより、頭頸部器官及び擬似経口摂取品１００を形成している粒子１つ１つを表示させる場合に比して、頭頸部器官及び擬似経口摂取品１００の表示形態を簡略化でき、開発者等に対して、頭頸部器官の運動や、擬似経口摂取品１００の挙動を容易に確認させることができる。
【0069】
（４）＜動的三次元頭頸部モデルの作製＞
このような動的三次元頭頸部モデル１０ｃは、特開２０２１－１１２５５７号公報で開示された作製手法に従って作製することができる。すなわち、医学的知見により理解されている頭頸部の構造、及び、ＣＴ（コンピュータ断層撮影：Computed Tomography）画像により大まかに読み取ることのできる口蓋１３、舌１２、気管１６の形態から、咽頭部１４と食道入口部１７ａの位置を推定し、舌１２や、口蓋１３、咽頭部１４、喉頭蓋１５ａ、喉頭１５等の頭頸部器官の構造を、ＣＧ（コンピュータグラフィックス：Computer Graphics）用ソフトウェア（Autodesk 3ds Max等）を用いてモデリングし、嚥下に関わる頭頸部器官を三次元（立体構造）で表した静的初期形状モデル（図示せず）を作製する。
このような動的三次元頭頸部モデル１０ｃは、特開２０２１－１１２５５７号公報で開示された作製手法に従って作製することができる。すなわち、医学的知見により理解されている頭頸部の構造、及び、ＣＴ（コンピュータ断層撮影：Computed Tomography）画像により大まかに読み取ることのできる口蓋１３、舌１２、気管１６の形態から、咽頭部１４と食道入口部１７ａの位置を推定し、舌１２や、口蓋１３、咽頭部１４、喉頭蓋１５ａ、喉頭１５等の頭頸部器官の構造を、ＣＧ（コンピュータグラフィックス：Computer Graphics）用ソフトウェア（Autodesk 3ds Max等）を用いてモデリングし、嚥下に関わる頭頸部器官を三次元（立体構造）で表した静的初期形状モデル（図示せず）を作製する。
【0070】
得られた静的初期形状モデルに対して、ＶＦ（嚥下造影検査：Videofluoroscopic examination of swallowing）による嚥下時の造影画像（正面及び側面図）を重ね合わせて、立体構造を修正し、被験者の安静時における頭頸部器官の立体的形状をＣＧによって描画した静的三次元頭頸部モデル（図示せず）を作製する。または、嚥下中の４次元ＣＴ（Computed Tomography）画像（４ＤＣＴ画像）をもとにして静的三次元頭頸部モデルを作成することもできる。このような静的三次元頭頸部モデルは、頭頸部モデリング部１０により作製される。次に、ＣＧで作製した静的三次元頭頸部モデルにおける頭頸部器官の表面を特定し、頭頸部器官の表面を境界として、頭頸部器官ごとにそれぞれの領域内に粒子を配置する。
【0071】
本実施形態における粒子は、三次元画像内で立体的形状を有した、三次元の球状粒子（本実施形態では単に粒子と称する）であり、例えば、乳児又は成人男性の平均的な大きさの頭頸部の原寸大モデルを三次元画像で作製する際には粒子の直径を0.1mm～3.0mm程度とすることが望ましく、そのうち、より好ましくは直径が0.6mm～1.5mm程度であることが望ましい。
【0072】
また、三次元画像内において作製した静的三次元頭頸部モデルにおいて、喉頭蓋の厚さ（例えば、成人では約3.0mm程度の厚さであり、乳児では約1.5mm程度の厚さである）方向に対して粒子が、少なくとも２個以上形成されることが望ましい。粒子の直径が小さすぎるとパーソナルコンピュータ２の計算処理負担が大きくなりすぎるため好ましくなく、一方、粒子の直径が大きすぎると、頭頸部器官について細かな運動を再現できないため、粒子の直径は上記の範囲とすることが望ましい。
【0073】
本実施形態においては、頭頸部器官を形成する粒子として、三次元の球状粒子を適用した場合について述べるが、本発明はこれに限らず、直方体形状の粒子、多角形状の粒子等その他種々の形状でなる粒子により頭頸部器官を形成するようにしてもよい。
【0074】
動的三次元頭頸部モデル１０ｃを作製する際は、ＣＴ画像及びＶＦ画像に基づいてＣＧにより作製した静的三次元頭頸部モデルの各頭頸部器官の表面が特定された後に、当該表面に囲まれた領域内に粒子同士が接するようにして粒子が隙間なく配置されることで、各頭頸部器官について粒子によるモデル化が行われる。
【0075】
また、動的三次元頭頸部モデル１０ｃで擬似経口摂取品を嚥下させる嚥下シミュレーションを行う場合には、擬似経口摂取品についても、ＣＧにより作製した擬似経口摂取品１００の表面が特定された後に、当該表面に囲まれた領域内に粒子同士が接するように粒子が隙間なく配置されることで、擬似経口摂取品１００について粒子によるモデル化が行われる。
【0076】
このようにして、静的三次元頭頸部モデルの各頭頸部器官及び擬似経口摂取品１００を、それぞれ粒子により作製する処理は、頭頸部モデリング部１０により行われる。
【0077】
ここで、本実施形態における動的三次元頭頸部モデル１０ｃでは、複数の粒子のうち、所定領域にある粒子を、擬似経口摂取品１００の嚥下時に頭頸部器官で強制的に移動する強制移動粒子として設定している。また、本実施形態における動的三次元頭頸部モデル１０ｃでは、複数の粒子のうち、強制移動粒子とした粒子以外で所定領域にある粒子を、筋線維方向に基づく収縮応力が与えられる筋粒子として設定している。
【0078】
なお、強制移動粒子及び筋粒子以外の粒子は、嚥下時に頭頸部器官で強制的に移動する位置（すなわち、嚥下時に三次元画像内で移動する座標）が規定されておらず、かつ、筋粒子のような筋線維方向への収縮応力についても規定されていない粒子である。このような強制移動粒子及び筋粒子以外の粒子は、動的三次元頭頸部モデル１０ｃで嚥下シミュレーションを行う際、従来の粒子法により移動位置等の解析が行われる。
【0079】
なお、強制移動粒子及び筋粒子以外の粒子の従来の嚥下シミュレーションの詳細については、文献「Kikuchi, T., Michiwaki, Y., Koshizuka, S., Kamiya, T., and Toyama Y., “Numerical simulation of interaction between organs and food bolus during swallowing and aspiration,” Computers in Biology and Medicine, 80, (2017), pp. 114‐123.」に開示されていることから、ここではその説明は省略し、以下、強制移動粒子と筋粒子とに着目して以下説明する。
【0080】
なお、次に説明する強制移動粒子に関し、粒子法を用いたシミュレーションについては、文献「Kikuchi, T., Michiwaki, Y., Kamiya, T. et al. Comp. Part. Mech. (2015) 2: 247. “Human swallowing simulation based on videofluorography images using Hamiltonian MPS method”」にも開示されている。
【0081】
動的三次元頭頸部モデル１０ｃは、頭頸部器官を形成する粒子の中に、強制移動粒子と筋粒子とを有している。動的三次元頭頸部モデル１０ｃは、頭頸部器官の一部の粒子を強制移動粒子として設定し、この強制移動粒子によって、嚥下時に主活動筋と考えられる筋の運動をモデル化している。この場合、動的三次元頭頸部モデル１０ｃでは、舌１２、口蓋１３、咽頭部１４、喉頭蓋１５ａ、喉頭１５及び食道１７等を形成している粒子の中から、嚥下時に、剛体的な強制変位を与える粒子を選定してこれを強制移動粒子１９として設定する。
【0082】
強制移動粒子１９は、実際に被験者が経口摂取品を嚥下する際における頭頸部器官の筋の運動が反映されるように、解剖学的知見や、医用画像の分析研究の知見から選定する。本実施形態では、所定の経口摂取品を被験者に嚥下させたときに得られたＶＦ画像や４ＤＣＴ画像において頭頸部器官をトレースし、動的三次元頭頸部モデル１０ｃ内で強制的に移動させる強制移動粒子１９を選定している。図４に示すように、動的三次元頭頸部モデル１０ｃには、舌強制移動部１９ａ、軟口蓋強制移動部１９ｂ、喉頭蓋強制移動部１９ｃ、甲状軟骨強制移動部１９ｄ、気管強制移動部１９ｅ、食道前壁強制移動部１９ｆを必要に応じて配置することができるが、必ずしもすべての強制移動部を用いる必要はなく、解析目的に応じて強制移動部を配置しないで実施することも可能である。
【0083】
また、所定の経口摂取品を被験者に嚥下させたときに得られたＶＦ画像や４ＤＣＴ画像に基づいて、嚥下開始から嚥下終了までの間、所定時間（例えば、０．１秒（以下、秒をＳとも記述する））ごとに各強制移動粒子１９が三次元画像内で移動する位置を決定し、各強制移動粒子１９について、嚥下時における時間と位置とを設定する。
【0084】
すなわち、動的三次元頭頸部モデル１０ｃにおいて、例えば、嚥下開始時である０．０Ｓのとき、三次元画像内のＸ軸、Ｙ軸及びＺ軸での座標が（０．０、０．２、－０．２）にある強制移動粒子１９について、嚥下開始時から０．１Ｓのとき座標（０．０、０．２、０．０）に移動し、０．２Ｓのとき座標（０．０、０．３、０．３）に移動することを設定する。なお、図３及び図４において黒丸で示した粒子は、本実施形態における強制移動粒子１９を示しており、例えば、舌１２、口蓋１３及び喉頭１５等の一部に設定されている。
【0085】
ここで、図９は、図４に示した動的三次元頭頸部モデル１０ｃにおいて、嚥下時に強制移動粒子１９が移動するときの軌跡の一部を、移動軌跡線２２ａ，２２ｂ，２２ｃ，２２ｄで表し、移動軌跡線２２ａ，２２ｂ，２２ｃ，２２ｄ（以下、これらをまとめて移動軌跡線２２とする）に従って、強制移動粒子１９を移動させたときの動的三次元頭頸部モデル１０ｃ１，１０ｃ２，１０ｃ３，１０ｃ４の状態変化を示した概略図である。
【0086】
例えば、２２ａは、舌１２に設定した強制移動粒子１９の移動軌跡線を示し、２２ｂは、口蓋１３の軟口蓋に設定した強制移動粒子１９の移動軌跡線を示し、２２ｃは、喉頭１５に設定した強制移動粒子１９の移動軌跡線を示し、２２ｄは、食道１７の管壁に設定した強制移動粒子１９の移動軌跡線を示す。
【0087】
図９では、約１２Ｓで所定の擬似経口摂取品を嚥下する動的三次元頭頸部モデル１０ｃを一例とし、嚥下開始時である０Ｓの動的三次元頭頸部モデル１０ｃ１と、嚥下開始から約７Ｓ後の動的三次元頭頸部モデル１０ｃ２と、嚥下開始から約９Ｓ後の動的三次元頭頸部モデル１０ｃ３と、嚥下開始から約１１Ｓ後の動的三次元頭頸部モデル１０ｃ４とを示す。
【0088】
このように、動的三次元頭頸部モデル１０ｃでは、頭頸部器官の所定位置に強制移動粒子１９を設け、各強制移動粒子１９が嚥下時に移動する位置と時間とを予め設定することで、嚥下時における頭頸部器官の基本的な運動（進行波的波動運動、回転運動、上下運動、前後運動、収縮運動等）を再現させている。なお、このような動的三次元頭頸部モデル１０ｃの強制移動粒子１９の運動は、器官運動設定部３０の強制運動設定部３１で設定する。
【0089】
動的三次元頭頸部モデル１０ｃにおいて嚥下時の咽頭部１４等の挙動を精度よく再現するためには、嚥下時に咽頭部１４等の壁面の長さが短縮する運動を再現することが望ましい。そこで、本実施形態における動的三次元頭頸部モデル１０ｃでは、咽頭部１４等の粒子に対して単に剛体的な強制変位を与えるだけでなく、咽頭部１４等の筋種ごとに三次元画像内で各筋粒子に筋線維方向を設定し、かつ筋線維方向に基づく最適な収縮応力を筋活動率により筋粒子に与え、動的三次元頭頸部モデル１０ｃにおける嚥下時の挙動を精度よく再現している。なお、動的三次元頭頸部モデル１０ｃにおいて筋線維方向に基づく収縮応力を与える粒子を筋粒子と称する。
【0090】
動的三次元頭頸部モデル１０ｃにおいて、例えば、舌１２を形成する粒子２０ａの中から筋粒子を設定する場合、図１０に示すように、解剖学的知見やＶＦ画像、４ＤＣＴ画像等に基づき、嚥下時に舌１２が伸縮する筋体領域ＥＲ１，ＥＲ２等を三次元画像内で特定し、各筋体領域ＥＲ１，ＥＲ２内に存在している粒子２０ａを探索する。例えば、筋体領域ＥＲ１内にある粒子２０ａを、舌１２の筋粒子とし、三次元画像の仮想空間内において、筋粒子ごとに筋線維方向を定義する。
【0091】
筋粒子に設定する筋線維方向の詳細については後述するが、解剖学的知見やＶＦ画像、４ＤＣＴ画像等に基づき、舌１２の筋体領域ＥＲ１内の空間内に、嚥下時に筋収縮が生じている方向を線分Ａとして複数設定し、筋粒子ごとに、近傍にある各線分Ａの方向の重み付け平均を筋線維方向としている。本実施形態においては、例えば、筋粒子から最も近い第１線分と、筋粒子に対して２番目に近い第２線分との２つの線分を特定し、筋粒子から第１線分までの距離と、筋粒子から第２線分までの距離とについて、筋粒子からの距離の近さに応じた重みを付けて第１線分の方向と第２線分の方向とを平均して筋線維方向を求めている。ただし、筋線維方向の求め方は、この手法である必要はなく、例えば、筋体領域ＥＲ１内に定義した全線分を用いて、放射基底関数（Radial Basis Function）補間を行うことでも、より滑らかに空間分布する、筋線維方向を得ることができる。
【0092】
また、動的三次元頭頸部モデル１０ｃにおいて、擬似経口摂取品を嚥下させたときに、頭頸部器官の収縮筋の筋種ごとに生じる、筋粒子の収縮応力の時間的変化を、筋活動率として設定し、筋活動率により嚥下時の収縮応力の大きさを設定している。
【0093】
ここでは、始めに、動的三次元頭頸部モデル１０ｃにおける咽頭部１４に着目し、咽頭部１４における筋種ごとに設定する筋線維方向について以下説明する。図１１の左図は、静的三次元頭頸部モデル１０ｂに、咽頭部１４の収縮筋が走行する方向Ａａ，Ａｂ，Ａｃ，Ａｄ，Ａｅ，Ａｆ，Ａｇ，Ａｈを示した筋線維モデル１０ｄの概略図である。図１１の右図は、左図の筋線維モデル１０ｄに示した咽頭収縮筋が走行する方向Ａａ，Ａｂ，Ａｃ，Ａｄ，Ａｅ，Ａｆ，Ａｇ，Ａｈを基に、動的三次元頭頸部モデル１０ｃにおける咽頭部１４において筋粒子１つ１つの筋線維方向を示した動的三次元頭頸部モデル１０ｅの概略図である。
【0094】
本実施形態では、筋線維モデル１０ｄ及び動的三次元頭頸部モデル１０ｅに示すように、解剖学的知見に基づき、上咽頭収縮筋舌咽頭部１４ａと、中咽頭収縮筋小角咽頭上部１４ｂと、中咽頭収縮筋小角咽頭下部１４ｃと、中咽頭収縮筋大角咽頭上部１４ｄと、中咽頭収縮筋大角咽頭下部１４ｅと、下咽頭収縮筋甲状咽頭上部１４ｆと、下咽頭収縮筋甲状咽頭下部１４ｇと、下咽頭収縮筋輪状咽頭部１４ｈとに、咽頭部１４を区分けしている。
【0095】
上咽頭収縮筋舌咽頭部１４ａと、中咽頭収縮筋小角咽頭上部１４ｂと、中咽頭収縮筋小角咽頭下部１４ｃと、中咽頭収縮筋大角咽頭上部１４ｄと、中咽頭収縮筋大角咽頭下部１４ｅと、下咽頭収縮筋甲状咽頭上部１４ｆと、下咽頭収縮筋甲状咽頭下部１４ｇと、下咽頭収縮筋輪状咽頭部１４ｈとについては、それぞれの領域内に筋粒子となる粒子が隙間なく配置され、粒子によるモデル化が行われている。
【0096】
そして、これら上咽頭収縮筋舌咽頭部１４ａと、中咽頭収縮筋小角咽頭上部１４ｂと、中咽頭収縮筋小角咽頭下部１４ｃと、中咽頭収縮筋大角咽頭上部１４ｄと、中咽頭収縮筋大角咽頭下部１４ｅと、下咽頭収縮筋甲状咽頭上部１４ｆと、下咽頭収縮筋甲状咽頭下部１４ｇと、下咽頭収縮筋輪状咽頭部１４ｈとについて、それぞれ収縮筋が走行する方向Ａａ，Ａｂ，Ａｃ，Ａｄ，Ａｅ，Ａｆ，Ａｇ，Ａｈは、解剖学的知見やＶＦ画像等に基づき、嚥下時に、各筋収縮筋の部位ごとにそれぞれ生じる細かな筋収縮方向を線分として設定した後、筋粒子ごとに、近傍にある各線分の方向の重み付け平均をして求めたものである。なお、図１１の左側に示した、収縮筋が走行する方向Ａａ，Ａｂ，Ａｃ，Ａｄ，Ａｅ，Ａｆ，Ａｇ，Ａｈは、説明の便宜上、筋線維方向のおおよその方向を示したものである。咽頭部１４の各筋粒子１つ１つは、このような収縮筋が走行する方向Ａａ，Ａｂ，Ａｃ，Ａｄ，Ａｅ，Ａｆ，Ａｇ，Ａｈに従って筋線維方向を定義している。
【0097】
上述した上咽頭収縮筋舌咽頭部１４ａと、中咽頭収縮筋小角咽頭上部１４ｂと、中咽頭収縮筋小角咽頭下部１４ｃと、中咽頭収縮筋大角咽頭上部１４ｄと、中咽頭収縮筋大角咽頭下部１４ｅと、下咽頭収縮筋甲状咽頭上部１４ｆと、下咽頭収縮筋甲状咽頭下部１４ｇと、下咽頭収縮筋輪状咽頭部１４ｈとに設けられた各筋粒子には、それぞれ嚥下時の時間経過に合わせて筋活動率が設定される。
【0098】
本実施形態における嚥下シミュレーション装置１は、開発者が、動的三次元頭頸部モデル１０ｃを構成する各頭頸部器官における筋活動率の値をそれぞれ変更しながら嚥下シミュレーションを繰り返し行い、嚥下シミュレーション時における動的三次元頭頸部モデル１０ｃの挙動を解析することで、嚥下シミュレーション時に動的三次元頭頸部モデル１０ｃが所望する動作を実行可能な筋活動率を決定する。
【0099】
動的三次元頭頸部モデル１０ｃにおける頭頸部器官の動作や、擬似経口摂取品の挙動を三次元画像内に表して嚥下シミュレーションを行なう際、動的三次元頭頸部モデル１０ｃにおける頭頸部器官の筋粒子を、ムーニー・リブリン（Mooney‐Rivlin）体として粒子法（例えば、ハミルトニアン粒子法：Hamiltonian MPS法）により運動解析部５０で解析する。
【0100】
なお、本実施形態に係る動的三次元頭頸部モデル１０ｃは、特開２０２１－１１２５５７号公報で開示された動的三次元頭頸部モデルを適用しており、例えば、内力及び外力を組み込んだときの支配方程式や、筋粒子ｉに加わる弾性力ｆｉ,ｅｌａｓｔｉｃ、筋粒子ｉに加わる接触力ｆｉ,ｃｏｎｔａｃｔ、筋粒子ｉに加わる擬似経口摂取品からの流体力ｆｉ,ｉｎｔｅｒａｃｔｉｏｎ、筋活動率の設定等については同じであるため、ここではその説明は省略する。
【0101】
（５）＜動的三次元頭頸部モデルを用いた筋シナジーの解析手法＞
次に、上述した動的三次元頭頸部モデル１０ｃを用いて、嚥下時における筋シナジーの解析手法について以下説明する。本実施形態では、（ｉ）筋活動行列Ｘを生成する筋活動行列生成工程と、（ｉｉ）動的三次元頭頸部モデル１０ｃの嚥下シミュレーションの結果から得られた各筋種の筋活動率を取得して非負値行列因子分解を実行する非負値行列因子分解演算処理工程と、（ｉｉｉ）時間分割を行う時間分割処理工程と、（ｉｖ）時分割筋活動行列、空間パターン行列及び時分割時間パターン行列を解析する行列解析工程と、（ｖ）時間分割処理工程で得られた時分割筋活動行列に対してデータ補間を行う補間処理工程と、（ｖｉ）各筋種の筋シナジーへの集約処理工程と、を実行する。以下、順番に説明する。
次に、上述した動的三次元頭頸部モデル１０ｃを用いて、嚥下時における筋シナジーの解析手法について以下説明する。本実施形態では、（ｉ）筋活動行列Ｘを生成する筋活動行列生成工程と、（ｉｉ）動的三次元頭頸部モデル１０ｃの嚥下シミュレーションの結果から得られた各筋種の筋活動率を取得して非負値行列因子分解を実行する非負値行列因子分解演算処理工程と、（ｉｉｉ）時間分割を行う時間分割処理工程と、（ｉｖ）時分割筋活動行列、空間パターン行列及び時分割時間パターン行列を解析する行列解析工程と、（ｖ）時間分割処理工程で得られた時分割筋活動行列に対してデータ補間を行う補間処理工程と、（ｖｉ）各筋種の筋シナジーへの集約処理工程と、を実行する。以下、順番に説明する。
【0102】
なお、上述した本実施形態においては、（ｉ）筋活動行列生成工程と、（ｉｉ）非負値行列因子分解演算処理工程と、（ｉｉｉ）時間分割処理工程と、（ｉｖ）行列解析工程と、（ｖ）補間処理工程と、（ｖｉ）集約処理工程と、を実行する場合について説明するが、本発明はこれに限らない。例えば、このうち、（ｉｉｉ）時間分割処理工程と、（ｉｖ）行列解析工程と、（ｖ）補間処理工程とは実行せずに、（ｉ）筋活動行列生成工程と、（ｉｉ）非負値行列因子分解演算処理工程と、（ｖｉ）集約処理工程と、だけを実行してもよいし、（ｉｉ）非負値行列因子分解演算処理工程は実行せずに、（ｉ）筋活動行列生成工程と、（ｉｉｉ）時間分割処理工程と、（ｉＶ）行列解析工程と、（Ｖ）補間処理工程と、（Ｖｉ）集約処理工程と、だけを実行してもよい。
【0103】
（５－１）＜筋活動行列生成工程＞
ここで、本実施形態に係る動的三次元頭頸部モデル１０ｃでは、特開２０２１－１１２５５７号公報で開示されている技術内容を基に、筋種ごとに筋活動率が設定されている（特開２０２１－１１２５５７号公報参照）。
ここで、本実施形態に係る動的三次元頭頸部モデル１０ｃでは、特開２０２１－１１２５５７号公報で開示されている技術内容を基に、筋種ごとに筋活動率が設定されている（特開２０２１－１１２５５７号公報参照）。
【0104】
筋シナジー解析処理部９３は、運動解析部５０によって行われる動的三次元頭頸部モデル１０ｃの嚥下シミュレーション時に、嚥下開始から嚥下終了までの各筋種の筋活動率をそれぞれ取得する。図１２Ａ及び図１２Ｂは、動的三次元頭頸部モデル１０ｃで嚥下シミュレーションを行った際に嚥下開始から嚥下終了までの各筋種の筋活動率を筋シナジー解析処理部９３により取得した結果を示す。各グラフの横軸は時間を示し、縦軸は筋活動率を示す。なお、本実施形態に係る筋シナジー解析処理部９３は、動的三次元頭頸部モデル１０ｃに設定した４０種の筋種の筋活動率を取得している。
【0105】
筋シナジー解析処理部９３は、運動解析部５０による嚥下時における頭頸部器官の運動の嚥下シミュレーション結果から、嚥下時における筋種ごとの筋活動率の経時変化を表した筋活動行列Ｘを生成する。筋活動行列Ｘは、筋種ごとの筋活動率と、嚥下開始から嚥下終了までのサンプリング時刻とを行列状に配置した数値データであり、例えば、図１３に示すように、筋活動率を取得した筋種を縦列に示し、筋活動率を取得した嚥下開始から嚥下終了までのサンプリング時刻を横行に示して各筋種の筋活動率の経時変化を表した筋活動行列Ｘを生成する。
【0106】
（５－２）＜非負値行列因子分解演算処理工程＞
次に、上述した筋活動行列Ｘに対して非負値行列因子分解を行い、筋活動行列Ｘを空間パターン行列Ｗと、時間パターン行列Ｈとに変換する非負値行列因子分解演算処理工程について説明する。
次に、上述した筋活動行列Ｘに対して非負値行列因子分解を行い、筋活動行列Ｘを空間パターン行列Ｗと、時間パターン行列Ｈとに変換する非負値行列因子分解演算処理工程について説明する。
【0107】
筋シナジー解析処理部９３には、複数の筋種が協調運動によって制御されるとした筋シナジー仮説に基づいて開発者が任意に選択した筋シナジーの個数を表す数値が入力部８１を介して入力され、当該筋シナジーの個数が設定される。
【0108】
筋シナジー解析処理部９３は、設定された筋シナジーの個数を用いて筋活動行列Ｘに対して非負値行列因子分解を実行することにより、当該筋活動行列Ｘを、筋シナジーごとに構成筋種ごとの筋活動度合いをそれぞれ表した空間パターン行列Ｗと、筋シナジーごとの活動度合いの時系列データ（経時変化のデータ）を表した時間パターン行列Ｈとに変換する。なお、本実施形態では、説明を簡単にするため、設定する筋シナジーの個数を９個として以下説明するが、筋シナジーの個数としてはこれに限らない。
【0109】
具体的には、筋活動行列Ｘに対して非負値行列因子分解を実行することにより得られる空間パターン行列Ｗは、例えば、４０種の筋種（例えば、オトガイ舌筋（１）、オトガイ舌筋（２）…）を縦列に示し、９個の筋シナジーを横行に示した行列状の数値データであり、９個の筋シナジーごとにそれぞれ筋種ごとの筋活動率を表している。
【0110】
図１４、図１５及び図１６は、空間パターン行列Ｗの横行に示す９個の筋シナジーＳｙｎ１～Ｓｙｎ９のデータの詳細構成について一例を示したものである。筋シナジーＳｙｎ１～Ｓｙｎ９は、４０種の筋種を縦軸に示し、筋シナジーＳｙｎ１～Ｓｙｎ９ごとに嚥下時における各筋種の相対的な筋活動率度合いを示したデータを横軸に示したデータである。
【0111】
筋シナジー解析処理部９３は、図１４、図１５及び図１６に示すようなデータ構成の空間パターン行列Ｗを開発者に対して提示することで、開発者に対して、筋シナジーＳｙｎ１～Ｓｙｎ９ごとに変化する筋種の筋活動率度合いの違いを確認させ、筋シナジーＳｙｎ１～Ｓｙｎ９ごとに協調運動によって制御される筋種を推測させることができる。
【0112】
筋活動行列Ｘに対して非負値行列因子分解を実行することにより得られる時間パターン行列Ｈは、例えば、筋シナジーを縦列にそれぞれ示し、筋活動率を測定した嚥下開始から嚥下終了までのサンプリング時刻を横行に示した行列状の数値データであり、各筋シナジーの活動度合いの経時変化を表している。
【0113】
図１７は、時間パターン行列Ｈにおいて筋シナジーＳｙｎ１～Ｓｙｎ９ごとにそれぞれ得られる、時間パターン９個を１つにまとめたグラフである。図１７では、嚥下開始から嚥下終了までの時間を横軸に示し、筋シナジーＳｙｎ１～Ｓｙｎ９の活動度合いを縦軸に示している。
【0114】
時間パターン行列Ｈでは、図１７に示すように、筋シナジーＳｙｎ１～Ｓｙｎ９ごとの時間パターンを得ることができるため、開発者に対して、筋シナジーＳｙｎ１～Ｓｙｎ９のうち、いずれの筋シナジーがどの時間帯で活動しているかを確認させることができる。
【0115】
（５－３）＜時間分割処理工程＞
次に、上述した筋活動行列Ｘを解析して、筋活動行列Ｘ及び時間パターン行列Ｈについて時間分割を行う時間分割処理工程について説明する。
次に、上述した筋活動行列Ｘを解析して、筋活動行列Ｘ及び時間パターン行列Ｈについて時間分割を行う時間分割処理工程について説明する。
【0116】
動的三次元頭頸部モデル１０ｃで嚥下動作を再現させる際には、一連の嚥下動作の中に、精度よく頭頸部器官の動作を再現させることが必要な特徴的な動作が存在する。その一方で、動的三次元頭頸部モデル１０ｃによる嚥下動作の再現の精度が低くても、嚥下時の頭頸部器官の動作全体を検証する際に影響が少ない頭頸部器官の動作も存在する。そのため、図１３に示した筋活動行列Ｘの中で、動的三次元頭頸部モデル１０ｃで嚥下動作を精度よく再現させるのに重要となる領域ＥＲ３が存在すると言える。
【0117】
ここで、単に筋活動行列Ｘに対して非負値行列因子分解を実行した場合、筋活動行列Ｘ全体が均一な精度で、空間パターン行列Ｗと時間パターン行列Ｈとに近似されるだけである。そこで、本実施形態では、嚥下時における各頭頸部器官の特徴的な動作を考慮しながら、筋活動行列Ｘについて時間分割したり、時間パターン行列Ｈに制約を設けたりすることが望ましい。ここでは、始めに筋活動行列Ｘを時間分割する時間分割処理工程について説明し、時間パターン行列Ｈに対して制約を設ける時間パターン解析処理工程は、後述する「（５―４）＜行列解析工程＞」で説明する。
【0118】
この場合、筋シナジー解析処理部９３は、嚥下時における動的三次元頭頸部モデル１０ｃの動作を再現した嚥下シミュレーションの結果に基づいて筋活動行列Ｘを得ると、当該筋活動行列Ｘを開発者に対して表示部４を介して提示する。これにより、嚥下シミュレーション装置１は、嚥下時における動的三次元頭頸部モデル１０ｃの動作を再現した嚥下シミュレーションの結果（動画像）と、得られた筋活動行列Ｘと、に基づいて、筋活動行列Ｘの中で、嚥下時における頭頸部器官の特徴的な動作を精度よく再現させるのに重要となる領域ＥＲ１を、後述する分割手法を基に開発者に対して特定させる。
【0119】
筋シナジー解析処理部９３は、開発者が筋活動行列Ｘ内で特定した重要となる領域ＥＲ１に基づいて、開発者に対して筋活動行列Ｘの中から分割する時間帯を指定させ、得られた分割指定情報に基づいて、例えば、図１８の１８Ａに示すように、筋活動行列Ｘを時系列に沿って複数の時間領域に分割し、筋活動行列Ｘから複数の時分割筋活動行列Ｘｎを生成する。本実施形態では、筋活動行列Ｘを時系列に沿って６つの時間領域に分割して６つの時分割筋活動行列Ｘ１～Ｘ６を生成している例を示す。
【0120】
本実施形態における分割手法としては、図５の５Ａに示したオトガイ舌筋（３）の活動がまだ小さい時間帯を時分割筋活動行列Ｘ１とし、オトガイ舌筋（３）のみが大きく活動し、かつ、オトガイ舌筋（４）は活動前の時間帯を時分割筋活動行列Ｘ２とし、オトガイ舌筋（４）も活動し、かつ、舌骨を挙上させる筋が活動前の時間帯を時分割筋活動行列Ｘ３とし、舌骨を挙上させ、徐々に咽頭を閉鎖していく時間帯を時分割筋活動行列Ｘ４とし、舌骨を挙上させる筋が活動終了する直前で、かつ、下咽頭収縮筋輪状咽頭部が活動する時間帯を時分割筋活動行列Ｘ５とし、咽頭を開大させ、舌骨も復位させる時間帯を時分割筋活動行列Ｘ６としている。
【0121】
なお、オトガイ舌筋（３）及びオトガイ舌筋（４）は、舌骨の位置制御に重要な部位であり、舌の形状に大きく影響する。また、下咽頭収縮筋輪状咽頭部は、食塊を食道入口部から食道へ送り込む役割を持つ部位である。なお、オトガイ舌筋や下咽頭収縮筋輪状咽頭部以外に、オトガイ舌骨筋、顎舌骨筋等の他の筋の活動に注目し、時間分割を行ってもよい。
【0122】
筋活動行列Ｘを複数の時分割筋活動行列Ｘｎに時間分割する場合は、図１９に示すように、動的三次元頭頸部モデル１０ｃを用いた嚥下シミュレーションにより得られたシミュレーション結果から舌骨の運動を解析したり、図２０の２０Ａから２０Ｄに示すように、動的三次元頭頸部モデル１０ｃを用いた嚥下シミュレーションにより得られたシミュレーション結果から測定した各筋種の筋活動率を解析したりし、解剖学や生理学などの基礎医学的知見やリハビリテーション学や耳鼻咽喉科学、歯科学、理学療法学、言語療法学等の臨床医学的知見等から頭頸部器官の特徴的動作を解析して決めてゆくことが望ましい。
【0123】
図１８の１８Ａに示すように、筋活動行列Ｘを時系列に沿って６つの時間領域に分割（時間分割）した場合には、当該筋活動行列Ｘに対して非負値行列因子分解を実行することで得られた時間パターン行列Ｈについても、筋活動行列Ｘを分割した同じ時間領域で、６つの時間領域に分割することができる。
【0124】
筋シナジー解析処理部９３は、筋活動行列Ｘを分割した６つの時間領域に基づいて、筋活動行列Ｘと同じ時間領域に、時間パターン行列Ｈを時系列に沿って複数の時間領域に分割し、時間パターン行列Ｈから複数の時分割時間パターン行列Ｈｎを生成する。本実施形態では、筋活動行列Ｘを６つの時間領域に分割した時分割筋活動行列Ｘ１～Ｘ６に合わせて、時間パターン行列Ｈを時系列に沿って同じ６つの時間領域に分割して６つの時分割時間パターン行列Ｈ１～Ｈ６を生成している例を示す。
【0125】
（５－４）＜行列解析工程＞
次に、上述した「（５－３）＜時間分割処理工程＞」で時間パターン行列Ｈを時間分割した時分割時間パターン行列Ｈ１～Ｈ６等を解析する行列解析工程について説明する。ここで、時分割時間パターン行列Ｈ１～Ｈ６は、例えば、図１７に示した時間パターン行列Ｈを時系列に沿って６つの時間領域に分割したデータとなる。
次に、上述した「（５－３）＜時間分割処理工程＞」で時間パターン行列Ｈを時間分割した時分割時間パターン行列Ｈ１～Ｈ６等を解析する行列解析工程について説明する。ここで、時分割時間パターン行列Ｈ１～Ｈ６は、例えば、図１７に示した時間パターン行列Ｈを時系列に沿って６つの時間領域に分割したデータとなる。
【0126】
筋シナジー解析処理部９３は、上述したように筋活動行列Ｘを分割した６つの時間領域に基づいて、時間パターン行列Ｈを同じ６つの時間領域に分割して６つの時分割時間パターン行列Ｈ１～Ｈ６を生成すると、時分割時間パターン行列Ｈ１～Ｈ６の時間領域ごとに、それぞれの筋シナジーＳｙｎ１～Ｓｙｎ９の活動度合いを特定する。
【0127】
例えば、図１９及び図２０の２０Ａに示すように、第１の時分割筋活動行列Ｘ１の時間領域は、０．８～１．１［Ｓ］であり、同じ０．８～１．１［Ｓ］の間の第１の時分割時間パターン行列Ｈ１の時間領域では、図１７に示すように、第１の筋シナジーＳｙｎ１の活動度合いだけが数値を有している。
【0128】
これにより、筋シナジー解析処理部９３は、図１８の１８Ｂに示すように、第１の時分割時間パターン行列Ｈ１で第１の筋シナジーＳｙｎ１の行を有効値「１」とし、その他の第２の筋シナジーＳｙｎ２から第９の筋シナジーＳｙｎ９までの行をそれぞれ無効値「０」として規定する。
【0129】
また、例えば、図１９及び図２０の２０Ａに示すように、第２の時分割筋活動行列Ｘ２の時間領域は、１．２［Ｓ］であり、同じ１．２［Ｓ］での第２の時分割時間パターン行列Ｈ２の時間領域では、図１７に示すように、第２の筋シナジーＳｙｎ２の活動度合いだけが数値を有している。
【0130】
これにより、筋シナジー解析処理部９３は、図１８の１８Ｂに示すように、第２の時分割時間パターン行列Ｈ２で第２の筋シナジーＳｙｎ２の行を有効値「１」とし、その他の第１の筋シナジーＳｙｎ１や第３の筋シナジーＳｙｎ３から第９の筋シナジーＳｙｎ９までの行をそれぞれ無効値「０」として規定する。このようにして、筋シナジー解析処理部９３は、第３の時分割時間パターン行列Ｈ３や、第５の時分割時間パターン行列Ｈ５、第６の時分割時間パターン行列Ｈ６についても同様に有効値、無効値を規定する。
【0131】
ここで、第４の時分割筋活動行列Ｘ４の時間領域は、例えば、図１９、図２０の２０Ｃ及び２０Ｄに示すように、１．４～１．９［Ｓ］であり、同じ１．４～１．９［Ｓ］の間の第４の時分割時間パターン行列Ｈ４の時間領域では、図１７に示すように、第４の筋シナジーＳｙｎ４から第７の筋シナジーＳｙｎ７の複数の活動度合いが数値を有している。
【0132】
ここで、図１７では、の第４の時分割時間パターン行列Ｈ４の時間領域である１．４～１．９［Ｓ］の間のうち、１．４～１．６［Ｓ］の間に第４の筋シナジーＳｙｎ４の活動度合いが数値を有し、１．５～１．８［Ｓ］の間に第５の筋シナジーＳｙｎ５の活動度合いが数値を有し、１．７～１．９［Ｓ］の間に第６の筋シナジーＳｙｎ６の活動度合いが数値を有し、１．９［Ｓ］に第７の筋シナジーＳｙｎ７の活動度合いが数値を有する。
【0133】
このため、筋シナジー解析処理部９３は、図１８の１８Ｂに示すように、第４の時分割時間パターン行列Ｈ４で、第４の筋シナジーＳｙｎ４から第７の筋シナジーＳｙｎ７の行でそれぞれ活動度合いが数値を有する時間帯を有効値「１」とし、その他の筋シナジー活動率が数値を有しない時間帯や、第１の筋シナジーＳｙｎ１から第３の筋シナジーＳｙｎ３までの行、第８の筋シナジーＳｙｎ８及び第９の筋シナジーＳｙｎ９の行をそれぞれ無効値「０」として規定する。その結果、第４の時分割時間パターン行列Ｈ４では、有効値「１」を規定した領域が対角的に段差状に規定されている。
【0134】
このようにして、筋シナジー解析処理部９３は、時分割時間パターン行列Ｈ１～Ｈ６の時間領域ごとに、それぞれの筋シナジーＳｙｎ１～Ｓｙｎ９の活動度合いを特定してゆき、時分割時間パターン行列Ｈ１～Ｈ６において、活動度合いを有する筋シナジーＳｙｎ１～Ｓｙｎ９の行とその時間帯を有効値「１」とし、当該活動度合いを有しない筋シナジーＳｙｎ１～Ｓｙｎ９の行を無効値「０」として規定してゆく。
【0135】
ここで、図２１の２１Ａに示すように、第１の時分割筋活動行列Ｘ１は、第１の空間パターン行列Ｗ１と、第１の時分割時間パターン行列Ｈ１に近似される。まず、第１の時分割筋活動行列Ｘ１を、行列の積の因数分解により、第１の空間パターン行列Ｗ１・第１の時分割時間パターン行列Ｈ１（すなわち、Ｘ１ ≒ Ｗ１・Ｈ１）に変換する。次に、第１の時分割時間パターン行列Ｈ１において初期有効値「１」を与え、非負値行列因子分解演算による反復計算を行い、第１の空間パターン行列Ｗ１及び第１の時分割時間パターン行列Ｈ１の値を更新する。時分割筋活動行列Ｘ１は、更新された第１の空間パターン行列Ｗ１と、更新された第１の時分割時間パターン行列Ｈ１の筋シナジーＳｙｎ１の行と、の積で表される。同様に、図２１の２１Ｃに示すように、第４の時分割筋活動行列Ｘ４は、非負値行列因子分解演算による反復計算により、第４～第７までの複数の空間パターン行列Ｗ４～Ｗ７と、第４の時分割時間パターン行列Ｈ４と、の積に近似により表され、第６の時分割筋活動行列Ｘ６は、非負値行列因子分解演算による反復計算により、第９の空間パターン行列Ｗ９と、第６の時分割時間パターン行列Ｈ６と、の積に近似により表される。
【0136】
図２１の２１Ｂに示すように、第２の時分割筋活動行列Ｘ２は、第２の空間パターン行列Ｗ２（＝ Ｘ２）と、第２の時分割時間パターン行列Ｈ２において有効値「１」を与えた筋シナジーＳｙｎ２の行と、の積で表される。同様に、第３の時分割筋活動行列Ｘ３は、第３の空間パターン行列Ｗ３（＝ Ｘ３）と、第３の時分割時間パターン行列Ｈ３において有効値「１」を与えた筋シナジーＳｙｎ３の行と、の積で表され、第５の時分割筋活動行列Ｘ５は、第８の空間パターン行列Ｗ８（＝ Ｘ５）と、第５の時分割時間パターン行列Ｈ５において有効値「１」を与えた筋シナジーＳｙｎ８の行と、の積で表される。また、第４の時分割筋活動行列Ｘ４は、対応する第４～第７までの複数の空間パターン行列Ｗ４～Ｗ７と、第４の時分割時間パターン行列Ｈ４の筋シナジーＳｙｎ４～Ｓｙｎ７の行と、の積で表される。
【0137】
（５－５）＜補間処理工程＞
次に、データ補間を行う補間処理工程について説明する。ここで、１．４～１．９［Ｓ］の時間領域内における第４の時分割筋活動行列Ｘ４及び第４の時分割時間パターン行列Ｈ４では、複数の筋シナジーＳｙｎ４～Ｓｙｎ７が含まれていることから、１．４～１．９［Ｓ］の時間領域の０．１［Ｓ］の時間刻みをさらに細かくし、動的三次元頭頸部モデル１０ｃによる嚥下シミュレーション時における各筋種の筋活動率を補間して、当該動的三次元頭頸部モデル１０ｃの動作をスムーズにすることが望ましい。
次に、データ補間を行う補間処理工程について説明する。ここで、１．４～１．９［Ｓ］の時間領域内における第４の時分割筋活動行列Ｘ４及び第４の時分割時間パターン行列Ｈ４では、複数の筋シナジーＳｙｎ４～Ｓｙｎ７が含まれていることから、１．４～１．９［Ｓ］の時間領域の０．１［Ｓ］の時間刻みをさらに細かくし、動的三次元頭頸部モデル１０ｃによる嚥下シミュレーション時における各筋種の筋活動率を補間して、当該動的三次元頭頸部モデル１０ｃの動作をスムーズにすることが望ましい。
【0138】
図２２は、第４の時分割筋活動行列Ｘ４の時間領域である１．４［Ｓ］から１．９［Ｓ］までの間のオトガイ舌筋種（１）～（６）や舌骨舌筋等の筋活動率について示すデータである。図２２に示すように、第４の時分割筋活動行列Ｘ４の時間領域は、０．１［Ｓ］の時間刻みであり、１．４［Ｓ］から１．５［Ｓ］までの間に筋活動率は規定されていない。そこで、筋シナジー解析処理部９３は、第４の時分割筋活動行列Ｘ４の時間領域では、１．４［Ｓ］から１．５［Ｓ］までの間の時間帯も線形補間によって筋活動率を補間する。
【0139】
その後、筋シナジー解析処理部９３は、時間刻み間の筋活動率を補間した第４の時分割筋活動行列Ｘ４に対して、筋シナジーＳｙｎ４～Ｓｙｎ７の分割数である４個の設定値で非負値行列因子分解を実行することで、第４の時分割筋活動行列Ｘ４を、新たに空間パターン行列Ｗ４～Ｗ７と第４の時分割時間パターン行列Ｈ４とに変換する。
【0140】
筋シナジー解析処理部９３は、以上のようにして第４の時分割筋活動行列Ｘ４の筋活動率を線形補間して新たに得らえた第４の時分割時間パターン行列Ｈ４を、残りの第１の時分割時間パターン行列Ｈ１から第３の時分割時間パターン行列Ｈ３と、第５の時分割時間パターン行列Ｈ５と、第６の時分割時間パターン行列Ｈ６とに統合して、図１７に示すように、筋シナジーＳｙｎ１～Ｓｙｎ９ごとに活動度合いの経時変化を表した時間パターンを新たに得る。
【0141】
なお、上述した実施形態においては、筋活動率を線形補間した第４の時分割筋活動行列Ｘ４に対して、非負値行列因子分解を実行することで、空間パターン行列Ｗ４～Ｗ７と第４の時分割時間パターン行列Ｈ４とを求めるようにした場合について説明したが、本発明はこれに限らない。例えば、第４の時分割筋活動行列Ｘ４の時間領域について筋活動率を線形補間した筋活動行列Ｘに対して非負値行列因子分解を実行することで、空間パターン行列Ｗと時間パターン行列Ｈとを求めるようにしてもよい。
【0142】
（５－６）＜集約処理工程＞
筋シナジー解析処理部９３は、図１７に示すように、筋シナジーＳｙｎ１～Ｓｙｎ９ごとに活動度合いの経時変化を表した時間パターンを得ると、４０種の筋種ごとに得られた各筋活動率を、筋種ごとに筋シナジーＳｙｎ１～Ｓｙｎ９の結果を反映した筋活動率へと集約する集約処理を実行する。
筋シナジー解析処理部９３は、図１７に示すように、筋シナジーＳｙｎ１～Ｓｙｎ９ごとに活動度合いの経時変化を表した時間パターンを得ると、４０種の筋種ごとに得られた各筋活動率を、筋種ごとに筋シナジーＳｙｎ１～Ｓｙｎ９の結果を反映した筋活動率へと集約する集約処理を実行する。
【0143】
具体的には、例えば、図１４に示すような第１の筋シナジーＳｙｎ１において、各筋種がそれぞれ第１の筋シナジーＳｙｎ１内で占める重みを算出する。第１の筋シナジーＳｙｎ１では、口蓋舌筋の占める重みが最も大きいことから、ここでは、口蓋舌筋に着目して以下説明する。
【0144】
筋シナジー解析処理部９３は、口蓋舌筋について第１の筋シナジーＳｙｎ１内で占める重みを算出すると、図１７に示すように、第１の筋シナジーＳｙｎ１の経時変化する活動度合いに対して、口蓋舌筋について第１の筋シナジーＳｙｎ１内で占める重みを乗算し、嚥下開始から嚥下終了までの口蓋舌筋の筋活動を第１の筋シナジーＳｙｎ１の活動度合いで表した第１データを得る。
【0145】
同様にして、筋シナジー解析処理部９３は、口蓋舌筋について第２の筋シナジーＳｙｎ２内で占める重みを算出し、図１７に示すように、第２の筋シナジーＳｙｎ２の経時変化する活動度合いに対して、口蓋舌筋について第２の筋シナジーＳｙｎ２内で占める重みを乗算し、嚥下開始から嚥下終了までの口蓋舌筋の筋活動を第２の筋シナジーＳｙｎ２の活動度合いで表した第２データを得る。
【0146】
このようにして、筋シナジー解析処理部９３は、筋シナジーＳｙｎ１～Ｓｙｎ９ごとに、嚥下開始から嚥下終了までの口蓋舌筋の筋活動率を各筋シナジーＳｙｎ１～Ｓｙｎ９の重みから算出した筋活動率でそれぞれ表した第１筋活動率から第９筋活動率を得る。得られた第１筋活動率から第９筋活動率は、筋シナジーの解析結果を反映した口蓋舌筋の筋活動率となる。
【0147】
筋シナジー解析処理部９３は、このようにして４０種の全ての筋種について第１筋活動率から第９筋活動率を算出し、それぞれの筋種ごとに筋シナジーＳｙｎ１～Ｓｙｎ９の結果を集約した筋活動率を算出する。このようにして、筋シナジー解析処理部９３は、４０種の筋種全てについて、筋シナジーＳｙｎ１～Ｓｙｎ９の解析結果を反映した筋活動率を得る。
【0148】
運動解析部５０は、４０種の筋種についてそれぞれ嚥下開始から嚥下終了まで経時変化する筋シナジーの解析結果を反映した筋活動率を動的三次元頭頸部モデル１０ｃに設定し、当該動的三次元頭頸部モデル１０ｃを用いて、嚥下時における頭頸部器官の運動を粒子法に基づいて三次元画像でシミュレーションする。
【0149】
ここで、筋シナジーの個数が適切な個数であれば、動的三次元頭頸部モデル１０ｃにおいて所望する嚥下動作を実現したシミュレーション結果が得られる。一方、筋シナジーの個数が少ない等の不適切な個数であれば、例えば、動的三次元頭頸部モデル１０ｃにおいて喉頭蓋が反転せずに喉頭に蓋をしない等の所望しない嚥下動作を実現したシミュレーション結果となる。
【0150】
なお、検証試験として、上述した実施形態に係る嚥下シミュレーション装置１を用いて、筋シナジーの個数を５個に設定し、上述した「（５）＜動的三次元頭頸部モデルを用いた筋シナジーの解析手法＞」に従って４０種の筋種全てについて、５個の筋シナジーＳｙｎ１～Ｓｙｎ５の筋シナジーの結果を反映した筋活動率を得、当該筋シナジーの解析結果を反映した筋活動率を設定した動的三次元頭頸部モデル１０ｃを用いて、嚥下時における頭頸部器官の運動を粒子法に基づいて三次元画像でシミュレーションしたところ、喉頭蓋が反転せずに気管に蓋をしない等の挙動となり、所望する動作を再現できなかった。そのため、動的三次元頭頸部モデル１０ｃで所望する嚥下動作を再現するためには、筋シナジーの個数を変更する必要があることが推察できた。
【0151】
（６）＜シミュレーション判定処理＞
シミュレーション判定部９４は、４０種の筋種の筋活動率を筋シナジーの解析結果を反映した筋活動率に変えた動的三次元頭頸部モデル１０ｃを用いて、運動解析部５０による嚥下シミュレーションを行った結果から、筋シナジーの個数の修正が必要か否かを判定する。具体的には、シミュレーション判定部９４は、動的三次元頭頸部モデル１０ｃの嚥下シミュレーションが最適な嚥下動作であるか否かを判定可能な機械学習モデル等を予め記憶しておく。シミュレーション判定部９４は、設定した筋シナジーの個数により動的三次元頭頸部モデル１０ｃの嚥下シミュレーションが最適な嚥下動作であるか否かを機械学習モデルにより判定する。
シミュレーション判定部９４は、４０種の筋種の筋活動率を筋シナジーの解析結果を反映した筋活動率に変えた動的三次元頭頸部モデル１０ｃを用いて、運動解析部５０による嚥下シミュレーションを行った結果から、筋シナジーの個数の修正が必要か否かを判定する。具体的には、シミュレーション判定部９４は、動的三次元頭頸部モデル１０ｃの嚥下シミュレーションが最適な嚥下動作であるか否かを判定可能な機械学習モデル等を予め記憶しておく。シミュレーション判定部９４は、設定した筋シナジーの個数により動的三次元頭頸部モデル１０ｃの嚥下シミュレーションが最適な嚥下動作であるか否かを機械学習モデルにより判定する。
【0152】
シミュレーション判定部９４は、最適な嚥下動作を再現していないとの判定結果を得たときには、筋シナジー解析処理部９３において設定する筋シナジーの個数の修正等を促し、修正した筋シナジーの個数によって、再び、筋活動行列Ｘに対して非負値行列因子分解を実行させるように促すことができる。
【0153】
このようにして、嚥下シミュレーション装置１では、筋種が協調運動する筋シナジーの個数や、嚥下時における筋シナジーごとの構成筋種の筋活動度合いを示した空間パターンと、筋シナジーごとの活動度合いの経時変化を示した時間パターンとを解析しながら、嚥下シミュレーションを行ってゆくことができる。
【0154】
（７）＜作用及び効果＞
以上の構成において、嚥下シミュレーション装置１は、複数の頭頸部器官が筋種ごとに三次元画像内で筋線維方向を特定した複数の筋粒子によりモデル化され、かつ、嚥下時における筋線維方向に基づく収縮応力が筋活動率により筋粒子に与えられる動的三次元頭頸部モデル１０ｃを、記憶部８３に記憶する。運動解析部５０は、動的三次元頭頸部モデル１０ｃを用い、嚥下時における頭頸部器官の運動を粒子法に基づいて三次元画像でシミュレーションする。
以上の構成において、嚥下シミュレーション装置１は、複数の頭頸部器官が筋種ごとに三次元画像内で筋線維方向を特定した複数の筋粒子によりモデル化され、かつ、嚥下時における筋線維方向に基づく収縮応力が筋活動率により筋粒子に与えられる動的三次元頭頸部モデル１０ｃを、記憶部８３に記憶する。運動解析部５０は、動的三次元頭頸部モデル１０ｃを用い、嚥下時における頭頸部器官の運動を粒子法に基づいて三次元画像でシミュレーションする。
【0155】
筋シナジー解析処理部９３は、運動解析部５０による嚥下時における頭頸部器官の嚥下運動のシミュレーション結果を基に、当該頭頸部器官の運動の経時変化に応じて変化する各筋種の筋活動率を取得する。筋シナジー解析処理部９３は、動的三次元頭頸部モデル１０ｃから取得した各筋種の筋活動率に基づいて、複数の筋種が協調運動によって制御されるとした筋シナジー仮説に基づいて定める筋シナジーの個数と、嚥下時における筋シナジーごとの構成筋種の筋活動度合いを示す空間パターンと、筋シナジーごとの活動度合いを経時変化で示した時間パターンと、を解析することができる。
【0156】
このように嚥下シミュレーション装置１は、筋種ごとに筋粒子に筋活動率が設定され当該筋活動率が変更可能な動的三次元頭頸部モデル１０ｃを用いて、嚥下時における頭頸部器官の運動を粒子法に基づいてシミュレーションし、当該シミュレーション時における筋活動率の解析結果に基づいて嚥下時における頭頸部器官の運動に関して筋シナジーの解析を行え得る。よって、嚥下シミュレーション装置１では、嚥下時における頭頸部器官の運動に関して、筋シナジーを含めたシミュレーションを行え得る。
【0157】
また、筋シナジー解析処理部９３は、運動解析部５０による嚥下時における頭頸部器官の運動のシミュレーション結果から、嚥下時における筋種ごとの筋活動率の経時変化を表した筋活動行列Ｘを生成する。そして、筋シナジー解析処理部９３は、設定された筋シナジーの個数を用いて筋活動行列Ｘに対して非負値行列因子分解を実行することにより、筋活動行列Ｘを空間パターン行列Ｗと時間パターン行列Ｈとに変換し、空間パターン及び時間パターンの解析結果として、これら空間パターン行列Ｗと時間パターン行列Ｈとを算出する。これにより、筋シナジー解析処理部９３は、筋活動行列Ｘを用いた演算処理によって空間パターン行列Ｗと時間パターン行列Ｈとを簡単に得られ、嚥下時における筋シナジーの空間パターンと時間パターンを解析することができる。
【0158】
さらに、筋シナジー解析処理部９３は、筋活動行列Ｘを時系列に沿って分割した時間領域に対応して、時間パターン行列Ｈも時系列に沿って当該時間領域に分割して複数の時分割時間パターン行列Ｈｎを生成する。そして、得られた複数の時分割時間パターン行列Ｈｎのうち、複数の筋シナジーＳｙｎ１～Ｓｙｎ９の活動度合いを含む時分割時間パターン行列Ｈ４の時間領域を特定して、特定した時間領域においてサンプリング時刻間の筋種の筋活動率を補間する。このように、嚥下シミュレーション時における各筋種の筋活動率を補間した分、動的三次元頭頸部モデル１０ｃの動作をスムーズにすることができる。
【0159】
さらに、筋シナジー解析処理部９３は、筋シナジーＳｙｎ１～Ｓｙｎ９ごとに得られた活動度合いを集約した筋シナジーの解析結果を反映した筋活動率を筋種ごとに生成し、筋種ごとに得られた筋シナジーの解析結果を反映した筋活動率を筋種ごとに設定する。これにより、運動解析部５０は、筋シナジー解析処理部９３で筋種ごとに得られた筋シナジーの解析結果を反映した筋活動率を筋種ごとに設定した動的三次元頭頸部モデル１０ｃにより、嚥下時における頭頸部器官の運動を粒子法に基づいて三次元画像でシミュレーションする。これにより、シミュレーション結果から筋シナジーの個数等を解析することができる。
【0160】
（８）＜他の実施形態＞
上述した実施形態においては、嚥下障害者又は健常者の頭頸部器官を再現した動的三次元頭頸部モデル１０ｃを三次元画像により形成する他にも、例えば、乳幼児又は高齢者等の頭頸部器官を再現した動的三次元頭頸部モデルを三次元画像により形成してもよい。また、上述した実施形態においては、粒子法として、本実施形態で採用したＭＰＳ法以外に、ＳＰＨ（Smoothed Particle Hydrodynamics）法や離散化要素法（Discrete Element Method）等を適用してもよい。
上述した実施形態においては、嚥下障害者又は健常者の頭頸部器官を再現した動的三次元頭頸部モデル１０ｃを三次元画像により形成する他にも、例えば、乳幼児又は高齢者等の頭頸部器官を再現した動的三次元頭頸部モデルを三次元画像により形成してもよい。また、上述した実施形態においては、粒子法として、本実施形態で採用したＭＰＳ法以外に、ＳＰＨ（Smoothed Particle Hydrodynamics）法や離散化要素法（Discrete Element Method）等を適用してもよい。
【0161】
また、上述した実施形態においては、動的三次元頭頸部モデルの頭頸部器官の運動（動作）や、擬似経口摂取品の挙動を、三次元画像内に表してシミュレーションする数値解析の方法として、粒子法を採用する場合について説明したが、本発明はこれに限らず、動的三次元頭頸部モデルの頭頸部器官や、擬似経口摂取品をメッシュとして扱う格子法を採用するようにしてもよい。この場合、粒子状の筋形成体である筋粒子に代わりに、メッシュ状の筋形成体（以下、単に筋メッシュと称する）を用いて動的三次元頭頸部モデル１０ｃを形成する。
【0162】
すなわち、動的三次元頭頸部モデル１０ｃとして、複数の頭頸部器官が筋種ごとに三次元画像内で筋線維方向を特定した、複数の筋メッシュによりモデル化され、かつ、嚥下時における筋線維方向に基づく収縮応力が筋活動率により筋メッシュに与えられる動的三次元頭頸部モデルを用いる。そして、運動解析部５０では、このような動的三次元頭頸部モデル（動的三次元頭頸部メッシュモデル）を用いて、格子法を用いた嚥下シミュレーションを行い、嚥下時における頸部器官の運動（動作）や、経口摂取品の嚥下時の挙動を、解析し、三次元画像内に表してもよい。なお、ここでの筋活動率は、嚥下時に筋線維方向に基づく収縮応力が与えられる筋メッシュに関し、筋種ごとに設定される、嚥下時における筋メッシュの収縮応力の時間的変化を示すものである。
【0163】
格子法としては、有限要素法（Finite Element Method）、有限体積法（Finite Volume Method）、有限差分法（Finite Difference Method）、又は境界要素法（Boundary Element Method）を適用することが望ましい。また、ＶＯＦ（Volume Of Fluid）法などの界面捕獲法を適用することも可能である。
【0164】
ここで、本実施形態に係る嚥下シミュレーション装置１にて使用する動的三次元頭頸部モデル１０ｃは、粒子法での説明に準じて作製することができる。すなわち、（４）＜動的三次元頭頸部モデルの作製＞と同様に静的三次元頭頸部モデルを作製後、頭頸部器官の表面を境界として、頭頸部器官ごとにそれぞれの領域内に格子を作製する。格子は、例えば、特開平４－１８１４８１号公報やＰａｌ ｅｔ．ａｌ．，Ｐｒｏｃ．Ｒ．Ｓｏｃ．Ｌｏｎｄ．Ｂ．，２７１，２５８７（２００４）等に開示された既知の手法に従って作製することができる。
【0165】
本実施形態における格子は、三次元画像内で立体的形状を有した、三次元の立方体格子又は直方体格子（本実施形態では単に格子と称する）であり、例えば、乳児又は成人男性の平均的な大きさの頭頸部の原寸大モデルを三次元画像で作製する際には格子の各辺の大きさを0.1mm～3.0mm程度とすることが望ましく、そのうち、より好ましくは各辺の大きさが0.6mm～1.5mm程度であることが望ましい。
【0166】
また、三次元画像内において作製した静的三次元頭頸部モデルにおいて、喉頭蓋の厚さ（例えば、成人では約3.0mm程度の厚さであり、乳児では約1.5mm程度の厚さである）方向に対して格子が、少なくとも２個以上形成されることが望ましい。格子が小さすぎるとパーソナルコンピュータ２の計算処理負担が大きくなりすぎるため好ましくなく、一方、格子が大きすぎると、頭頸部器官について細かな運動を再現できないため、格子の大きさは上記の範囲とすることが望ましい。
【0167】
本実施形態においては、頭頸部器官を形成する格子として、三次元の立方体又は直方体を適用した場合について述べるが、本発明はこれに限らず、三角形状の格子、多角形状の格子等その他種々の形状でなる格子により頭頸部器官を形成するようにしてもよい。
【0168】
このように、上述した実施形態の筋粒子を筋メッシュに変えてモデル化した動的三次元頭頸部モデルを用いる嚥下シミュレーション装置でも、上述した実施形態と同様に、筋シナジー解析処理部９３によって、当該動的三次元頭頸部モデルから取得した筋種ごとの筋メッシュの筋活動率に基づいて、複数の筋種が協調運動によって制御されるとした筋シナジーの個数と、嚥下時における筋シナジーごとの構成筋種の筋活動度合いを示す空間パターンと、筋シナジーごとの活動度合いを経時変化で示した時間パターンと、を解析することができる。
【0169】
このように嚥下シミュレーション装置でも、同様に、筋種ごとに筋メッシュに筋活動率が設定され当該筋活動率が変更可能な動的三次元頭頸部モデルを用いて、嚥下時における頭頸部器官の運動を格子法に基づいてシミュレーションし、当該シミュレーション時における筋活動率の解析結果に基づいて嚥下時における頭頸部器官の運動に関して筋シナジーの解析を行え得る。よって、嚥下シミュレーション装置では、嚥下時における頭頸部器官の運動に関して、筋シナジーを含めたシミュレーションを行え得る。
【0170】
上述した実施形態においては、シミュレーション装置及びシミュレーション方法として、擬似経口摂取品の嚥下時における頭頸部器官の運動を粒子法又は格子法に基づいて三次元画像でシミュレーションする嚥下シミュレーション装置１及び嚥下シミュレーション方法について説明したが、本発明はこれに限らず、上述したように、擬似経口摂取品の咀嚼時における頭頸部器官の運動を粒子法又は格子法に基づいて三次元画像でシミュレーションする咀嚼シミュレーション装置及び咀嚼シミュレーション方法に適用してもよい。このような、咀嚼シミュレーション装置でも、上述した実施形態と同様に、咀嚼時における頭頸部器官の運動に関して、筋シナジーを含めたシミュレーションを行え得る。
【符号の説明】
【0171】
１ 嚥下シミュレーション装置（シミュレーション装置）
１０ｃ 動的三次元頭頸部モデル
５０ 運動解析部
８３ 記憶部
９３ 筋シナジー解析処理部
１０ｃ 動的三次元頭頸部モデル
５０ 運動解析部
８３ 記憶部
９３ 筋シナジー解析処理部
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12A】
【図12B】
【図13】
【図14】
【図15】
【図16】
【図17】
【図18】
【図19】
【図20】
【図21】
【図22】