JP6835014B2 - 身体内部情報推定方法、コンピュータプログラム、それを記憶した記憶媒体、および、身体内部情報推定装置 - Google Patents

Description

本発明は、身体内部情報推定方法、コンピュータプログラム、それを記憶した記憶媒体、および、身体内部情報推定装置に関する。

従来から、人体モデルを表すデータを用いて、人体内部に生じる力や人体内部の形状の変化を推定する技術が知られている（特許文献１〜３参照）。例えば、特許文献１には、骨格を表す剛体リンクを表すデータと、筋・腱・靱帯を表すデータにより定義された身体モデルに対して、身体代表部位の運動の加速度データを与え、逆運動学計算により筋・腱・靱帯の発生力を計算する技術が開示されている。特許文献２には、予め作成された人体内部組織の構造とその形状データを有する数値人体モデルの体表および内部組織の形状を変形させて、医療用のＣＴ／ＭＲＩ画像から得られた個人の体表や内部組織に近似させることにより、身体内部構造の３次元的な形状変化や配置変化を算出する技術が開示されている。特許文献３には、ＮＵＲＢＳ（Non-Uniform Rational B-Spline）による人体モデルの形状変形手法を用いて、人体の体表面の変形や、体表や内部に発生する力を算出する技術が開示されている。

特開２００３−３３９６７３号公報 特開２０１３−８９１２３号公報 米国特許出願公開第２０１０／０１５６９３５号公報

運動器系・循環器系の統合的数理モデルの構築：バイオメカニズム １８ システムとしての生体 pp.57-68,2006 三相理論を用いた心筋細胞の電気化学・力学連成シミュレーションの大規模化に関する検討：日本機械学会論文集（Ａ編）,79巻,803号,pp.934-949,2013

しかしながら、上記先行技術によっても、身体モデルによる身体運動から身体内部の情報を推定する技術については、なお、改善の余地があった。例えば、上記先行技術によっても、身体運動情報から血圧や血管・神経の伝達量などの身体内部情報を推定することは困難であった。

本発明は、上述した課題を解決するためになされたものであり、身体モデルを用いて、身体運動情報から身体内部情報を推定する技術の向上を図ることを目的とする。

本発明は、上述の課題の少なくとも一部を解決するためになされたものであり、以下の形態として実現することが可能である。

（１）本発明の一形態によれば、骨格を表す剛体リンクモデルのデータと、前記剛体リンクモデルに関連づけられた、骨、筋肉、血管、および、心臓を含む身体内部組織のモデルによって構成される身体モデルのデータと、を用いて、血圧や筋力を含む身体内部情報を推定する身体内部情報推定方法が提供される。この身体内部情報推定方法は、前記剛体リンクモデルに運動を表現させるための身体運動情報を入力する入力工程と、前記身体運動情報に基づく前記剛体リンクモデルの運動に追従して前記身体内部組織の形状変化を算出する変形計算工程と、数理モデルと数値解析の少なくとも一方を用いて、前記身体内部組織の形状変化から前記身体内部情報を推定する身体内部情報推定工程と、推定された前記身体内部情報を出力する出力工程と、を備える。

この構成によれば、剛体リンクモデルの運動に追従して身体内部組織の形状変化を算出する変形計算工程と、身体内部組織の形状変化から身体内部情報を推定する身体内部情報推定工程とを備えているため、身体運動情報から身体内部情報を推定することができる。

（２）上記形態の身体内部情報推定方法において、前記変形計算工程では、前記剛体リンクモデルの運動に追従して前記身体モデルをＦＦＤによって変形させることにより前記身体内部組織の形状変化を算出してもよい。この構成によれば、ＦＦＤ（Free Form Deformation）によって、身体モデルを変形させることにより、簡易に身体内部組織の形状変化を算出することができる。

（３）上記形態の身体内部情報推定方法において、前記身体内部情報推定工程では、数理モデルを用いて、前記身体内部情報のうちの血圧と筋力を算出するマクロ計算工程と、数値解析によって、前記身体内部情報のうちの血管血流量と神経電流量を算出するミクロ計算工程と、の少なくとも一方の工程が含まれていてもよい。この構成によれば、マクロ計算工程によって低計算コストで身体内部情報を算出できるとともに、ミクロ計算工程によって精確な身体内部情報を算出できる。

（４）上記形態の身体内部情報推定方法において、前記身体内部情報推定工程には、前記ミクロ計算工程において数値解析によって算出された前記身体内部情報を、前記マクロ計算工程において数理モデルを用いて前記身体内部情報を算出する際に利用する統合工程が含まれていてもよい。この構成によれば、数理モデルを用いて算出される身体内部情報の精度を向上させることができる。

（５）上記形態の身体内部情報推定方法において、 前記身体内部情報推定工程には、前記マクロ計算工程において数値モデルを用いて算出された前記身体内部情報を、前記ミクロ計算工程において数値解析によって前記身体内部情報を算出する際の境界条件として利用する統合工程が含まれていてもよい。この構成によれば、数値解析によって算出される身体内部情報の計算コストを低減させることができる。

（６）上記形態の身体内部情報推定方法において、前記身体内部情報推定工程には、前記変形計算工程で計算された特定の筋肉の形状変化から、前記筋肉の長さの変化と、前記筋肉の収縮速度を算出し、算出した前記筋肉の長さの変化と、前記筋肉の収縮速度から数理モデルを用いて前記特定の筋肉の筋力を算出する工程が含まれていてもよい。この構成によれば、筋肉の長さの変化、および、筋肉の収縮速度と、筋力との関係が示された数理モデルを用いて容易に筋力を算出することができる。

（７）上記形態の身体内部情報推定方法において、前記身体内部情報推定工程には、前記変形計算工程で計算された心臓の形状変化から、左心室の容積変化を算出し、算出した前記左心室の容積変化から数理モデルを用いて血圧を推定する工程が含まれていてもよい。この構成によれば、左心室の容積変化と血圧との関係が示された数理モデルを用いて容易に血圧を算出することができる。

（８）上記形態の身体内部情報推定方法において、前記身体内部情報推定工程には、前記身体内部組織のうち、血管と脂肪とが一つの混合体として構成された混合体モデルを用いて、有限要素解析により血管血流量を推定する工程が含まれていてもよい。この構成によれば、血管と脂肪を含む混合体において、時間ステップごとの平衡方程式の解を求めることにより血管血流量を算出することができる。

なお、本発明は、種々の態様で実現することが可能であり、例えば、身体内部情報推定装置、身体内部情報計算システム、身体内部情報推定方法、身体内部情報の推定をコンピュータに実行させるコンピュータプログラム、これを記憶した記憶媒体などの形態で実現することができる。

第１実施形態の身体内部情報推定装置の構成を例示した説明図である。 剛体リンクモデルを例示した説明図である。 人体モデルの全身を例示した説明図である。 人体モデルの腕部分を例示した説明図である。 人体モデルの胸部付近を例示した説明図である。 筋骨格モデルを例示した説明図である。 人体モデルと筋骨格モデルとの対応関係を示した説明図である。 ＦＦＤによる腕モデルの曲げを示した説明図である。 ＦＦＤによる心臓モデルのねじれ変形を示した説明図である。 身体内部情報計算部の構成を示した説明図である。 人の腕の断面と混合法用腕モデルを例示した説明図である。 第２実施形態の身体内部情報計算部の構成を例示した説明図である。 羽状筋を例示した説明図である。

＜第１実施形態＞
図１は、第１実施形態の身体内部情報推定装置１０の構成を例示した説明図である。身体内部情報推定装置１０は、身体モデルを用いたシミュレーションによって、身体運動情報から身体内の力学的・電気化学的情報（身体内部情報）を推定するコンピュータシステムである。身体内部情報推定装置１０は、演算処理装置、記憶装置、および、入出力装置を含んで構成され、記憶装置に記憶されたプログラムを演算処理装置が実行することにより後述の２つの計算部（身体運動追従変形計算部１２および身体内部情報計算部１３）を構成する。記憶装置には、このプログラムのほか、後述する種々の身体モデル（剛体リンクモデル、筋骨格モデル、人体モデル、混合法用モデル等）を表すデータが記憶されている。身体内部情報推定装置１０は、これら演算処理装置、記憶装置、および、入出力装置によって、身体運動情報入力部１１と、身体運動追従変形計算部１２と、身体内部情報計算部１３と、身体内部情報出力部１４と、が構成されている。

身体運動情報入力部１１は、入力装置によって構成され身体運動情報が入力される。身体運動情報とは、所定期間における人の運動軌跡を特定可能な情報であり、例えば、人の各部位の位置座標や速度、加速度等の時刻歴データや、筋活性度の時刻歴データ等が含まれる。この身体運動情報を剛体リンクモデルや筋骨格モデルに与えることによって、これらのモデルに身体運動を表現させることができる。身体運動情報、剛体リンクモデル、および、筋骨格モデルについては後述する。

身体運動追従変形計算部１２は、身体運動情報入力部１１に入力された身体運動情報を用いて身体運動を表現した際に、運動に追従して変形する身体内外部組織（皮膚、筋肉、血管、神経等）の形状変化を算出する。ここでは、身体運動情報を用いて身体運動を表現した際に、身体運動中の剛体リンクモデルに追従して、剛体リンクモデルの剛体部（例えば、上腕骨など）に関連づけられる皮膚、筋肉、血管、神経等の身体内外部組織を動かし、これらの形状を変化させる。形状を変化させる手法としては、例えば、コンピュータグラフィックス分野で用いられるＦＦＤ（Free Form Deformation）やＮＵＲＢＳ（Non-Uniform Rational B-Spline）を例示することができる。ＦＦＤによる身体内外部組織の形状変化の一例については後述する。

身体内部情報計算部１３は、身体運動追従変形計算部１２によって算出された身体内外部組織の変形時の形状変化に基づいて身体内部情報を計算する。身体内部情報とは、身体内部の力学的情報または電気化学的情報であり、例えば、筋力、血流量、血圧、神経電流量、応力、ひずみ等を例示することができる。これらの計算方法については後述する。

身体内部情報出力部１４は、出力装置によって構成され、身体内部情報計算部１３から得られた身体内部情報を出力する。

図２は、記憶装置に記憶されている剛体リンクモデル２０を例示した説明図である。剛体リンクモデル２０は、人体の骨を剛体のリンク２１とし、人体の関節部をジョイント２２としてモデル化したものである。身体運動追従変形計算部１２（図１）は、例えば、剛体リンクモデル２０に対して、モーションキャプチャー等によって計測された各関節や特徴点位置の運動データ（位置、速度、加速度）を関節角度や特徴点位置の時刻歴データとして与えることにより身体運動を表現することができる。

図３〜図５を用いて、記憶装置に記憶されている人体モデル３０について説明する。図３は、人体モデル３０の全身を例示した説明図である。図４は、人体モデル３０の腕部分を例示した説明図である。図５は、人体モデル３０の胸部付近を例示した説明図である。図３〜図５の全体で表される人体モデル３０は、骨３１、関節３２、皮膚３３、筋肉３４、血管３５、神経３６、心臓３７を含んで構成されている。図３には、人体モデル３０のうちの骨３１、関節３２、血管３５が示されている。図４には、人体モデル３０の腕部（腕モデル）３０ａにおける筋肉３４、血管３５、神経３６が示されている。図５（ａ）には、人体モデル３０の胸部における骨３１、皮膚３３、心臓３７が示されている。図５（ｂ）には、人体モデル３０の心臓３７のうちの左心室３７ａと、右心室３７ｂと、冠動脈３５ａと、大動脈３５ｂとが示されている。

ここでは、人体モデル３０の骨３１、関節３２は、剛体リンクモデル２０（図２）のリンク２１、ジョイント２２とそれぞれ関連づけられている。また、人体モデル３０の皮膚３３、筋肉３４、血管３５、神経３６、心臓３７は、剛体リンクモデル２０のリンク２１と関連づけられている。身体運動追従変形計算部１２（図１）は、例えば、剛体リンクモデル２０により身体運動を表現したときに、剛体リンクモデル２０のリンク２１に関連づけられている人体モデル３０の皮膚３３、筋肉３４、血管３５、神経３６、心臓３７を後述するＦＦＤ等によって変形させることができる。

図６〜図７を用いて、記憶装置に記憶されている筋骨格モデル４０について説明する。図６は、筋骨格モデル４０を例示した説明図であり、図７は、人体モデル３０と筋骨格モデル４０との対応関係を示した説明図である。図７では、人体モデル３０の一部である腕モデル３０ａと、筋骨格モデル４０の一部を構成する筋肉モデル４１が示されている。図６の筋骨格モデル４０は、図３の人体モデル３０から作成することができる。具体的には、図７に示すように、人体モデル３０（ここでは、腕モデル３０ａ）に含まれる各筋肉３４の３次元的な筋形状（筋ボリューム）を、両端の骨との付着点を結ぶ直線または曲線としてモデル化した筋肉モデル４１を作成し、これを図２の剛体リンクモデル２０に追加することによって筋骨格モデル４０を作成できる。このように、筋骨格モデル４０は、剛体リンクモデル２０と対応づけられており、この筋骨格モデル４０の筋肉モデル４１に、筋電計測などにより推定された筋活性度や筋制御アルゴリズムにより算出された筋活性度の時刻歴データを与えることにより、身体運動を表現することができる。

図８は、ＦＦＤによる腕モデル３０ａの曲げを示した説明図である。人体モデル３０の一部分である腕モデル３０ａの全体をＦＦＤ格子７１で囲み、肘関節を中心として前腕部に回転運動を与えることにより腕を屈曲させることができる。この場合、筋肉３４だけでなく内部に含まれる血管３５や神経３６も屈曲角度に応じて変形する。身体運動追従変形計算部１２（図１）は、剛体リンクモデル２０（図２）により身体運動を表現した際に、各時間ステップにおいてリンク２１やジョイント２２の位置や姿勢の変化に対応させてＦＦＤ格子を変形させる。これにより、身体内外部組織の形状変化をシミュレーションすることができる。このシミュレーションによって、例えば、運動時や運動前後におけける筋肉の長さの変化や筋肉の断面積の変化を算出することができる。

図９は、ＦＦＤによる心臓モデル３７のねじれ変形を示した説明図である。人体モデル３０の一部分である心臓モデル３７は、内部の左心室３７ａも表現されている。左心室３７ａを含む心臓モデル３７全体をＦＦＤ格子で囲み、Ｚ軸方向の収縮と、ＸＹ軸平面の収縮と、Ｚ軸回りのねじりと、の３つの変形に分けてＦＦＤ格子に与えることによって、運動軌跡として心拍一周期における左心室容積の時系列変化を与えることができる。身体運動追従変形計算部１２（図１）は、ＦＦＤを用いた心臓モデル３７の形状変化をシミュレーションすることにより、左心室の容積の変化を算出することができる。

図１０は、身体内部情報計算部１３の構成を示した説明図である。ここでは、身体内部情報計算部１３による身体内部情報の算出方法について説明する。身体内部情報計算部１３は、マクロ計算部１３１と、ミクロ計算部１３２とを備えている。

マクロ計算部１３１は、身体運動追従変形計算部１２によって算出された身体内外部組織の変形時の形状変化から、身体内の部位ごとの力学的・電気化学的情報をマクロ計算によって算出する。マクロ計算部１３１は、運動器系・循環器系の統合的数理モデルを使って、例えば、筋力、血流量、血圧、骨への力等を算出することができる。このマクロ計算は、身体運動や心臓の動きなどを低計算コストでシミュレーションできるため、学習の繰り返し計算に利用することができる。

ミクロ計算部１３２は、身体内外部組織の変形時の形状変化から、身体の組織ごとの力学・電気化学的情報をミクロ計算によって算出する。ミクロ計算部１３２は、数値解析（有限要素解析等）によって、例えば、血管血流量、神経電流量、筋・血管・神経に作用する応力・歪み、筋断面積変化、神経伝達物質の濃度等を算出することができる。このミクロ計算は、血管や神経に作用する力、各組織の詳細な変形形状、血流や神経伝達などの情報を精確に得るために利用することができる。

ここでは、マクロ計算の例として、運動器系・循環器系の統合的数理モデルを使った筋力、血流量、血圧の算出例について説明する。また、ミクロ計算の例として、混合法による有限要素解析を使った血管血流量、神経電流量の算出例について説明する。なお、マクロ計算、ミクロ計算の両者は相互に補間しながら利用してもよい。これについては第２実施形態で例示する。

＜筋力の算出例＞
身体運動追従変形計算部１２によって算出された身体内外部組織の変形時の形状変化から、身体内部情報計算部１３のマクロ計算部１３１が数理モデルを用いて筋力を算出する方法について説明する。マクロ計算部１３１は、第Ｍ筋の筋張力Ｆ^Mを下記の式（１）によって算出することができる。

式（１）において、ＰＣＳＡ^Mは筋断面積であり、ここでは、筋ごとに一定値と仮定する。α^Mは筋活性度であり０〜１００％の値をとる。σ^M _maxは単位筋断面積あたりの最大発揮筋力を表す係数である。Ｆ_CE,L（Ｌ^M）は能動的な筋力と筋の長さとの関係式である。Ｌ^Mは筋の自然長で正規化された筋長である。Ｆ_CE,V（Ｌ^M,Ｌ^*M）は能動的な筋力と速度との関係式である。Ｌ^*Mは最大収縮速度で正規化された筋収縮速度である。Ｆ_PE（Ｌ^M）は受動的な筋力と筋の長さとの関係式である。

式（１）に示すように、筋力Ｆ^Mは、筋断面積ＰＣＳＡ^Mと、単位面積あたりの力である筋応力との積で表現できる。式（１）の｛｝内（かっこ内）の第１項は、能動的な筋応力を示し、第２項は、受動的な筋応力を示している。各筋応力はその最大値σ^M _max（筋種毎に一定と仮定する、人の場合：５ｇ／ｃｍ²）と、筋長さおよび筋収縮速度に関係する変数との積で表現される。Ｆ_CE,L（Ｌ^M）は筋の長さ変化に依存する。Ｆ_CE,V（Ｌ^M,Ｌ^*M）は、筋の収縮速度に依存する。Ｆ_PE（Ｌ^M）は筋の長さ変化に依存する。α^Mは別途、筋制御アルゴリズムや筋電計測データから推定できる。

このように、筋力Ｆ^Mは、筋の長さ変化と、筋の収縮速度とを変数とする関数によって算出できる。時間ステップ毎の筋の長さ変化や筋の収縮速度は、身体運動追従変形計算部１２において筋骨格モデル４０やＦＦＤ変形計算から得られた腕の屈曲後の筋長さを測定することにより得ることができる。ＦＦＤ変形計算で腕を屈曲させる方法は図８で説明した通りである。筋骨格モデル４０で腕を屈曲させる場合、筋肉モデルに筋活性度の時刻歴データを与えて筋肉を収縮させるか、骨格部分である剛体リンクモデルに運動軌跡データを与えることによって屈曲させて、筋の長さを変化させることができる。以上のことから、身体内部情報計算部１３は、身体運動追従変形計算部１２によって算出された身体の変形からから筋力Ｆ^Mを算出することができる。

＜血流量、血圧の算出例＞
身体運動追従変形計算部１２によって算出された身体内外部組織の変形時の形状変化から、身体内部情報計算部１３のマクロ計算部１３１が数理モデルを用いて血流量および血圧を算出する方法について説明する。マクロ計算部１３１は、動脈の血圧Ｐ（ｔ）を下記の式（２）により算出することができる。また、心拍出量Ｑ（ｔ）を下記の式（３）により算出することができる。

式（２）において、Ｒは末梢血管抵抗である。θは大動脈のコンプライアンスであり、ここでは定数とする。Ωは１分あたりの心拍出量（ｍｌ）である。Ｔは心周期である。時刻ｔにおける心拍出量Ｑ（ｔ）（ｍｌ／ｓ）は左心室容積の時系列変化ｆ（ｔ／Ｔ）から得ることができる。すなわち、Ｑ（ｔ）は、身体運動追従変形計算部１２によるＦＦＤを用いた心臓モデル３７（図９）変形時の形状変化にともなう左心室３７ａの容積変化から算出することができる。以上のことから、身体内部情報計算部１３は、身体運動追従変形計算部１２によって算出された心臓の変形からから血圧Ｐ（ｔ）および心拍出量Ｑ（ｔ）を算出することができる。

＜血管血流量、神経電流量の算出例＞
身体運動追従変形計算部１２によって算出された身体内外部組織の変形時の形状変化から、身体内部情報計算部１３のミクロ計算部１３２が数値解析によって血管血流量および神経電流量を算出する方法について説明する。ここでは、図１１に示す混合法用モデルを用いた混合法による有限要素解析について説明する。

図１１は、人の腕の断面Ｓａと混合法用腕モデル５０を例示した説明図である。図１１（ａ）に示すように、腕の断面Ｓａには、皮膚５３、筋肉５４、骨５１がそれぞれ配置され、皮膚５３、筋肉５４、骨５１の隙間部分に、脂肪５２と神経５６と血管５５とを含む混合体５７が配置されている。筋肉５４は、上腕筋５４１と、上腕二頭筋５４２（長頭５４２ａおよび短頭５４２ｂ）と、上腕三頭筋５４３（長頭５４３ａ、外側頭５４３ｂ、および、内側頭５４３ｃ）と、を含んでいる。骨５１は、上腕骨である。図１１（ｂ）に示す混合法用腕モデル５０では、皮膚５３、筋肉５４、骨５１は、それぞれ有限要素法でメッシュ分割したモデルとして構成され、脂肪５２と神経５６と血管５５は、一つの混合体５７として有限要素法でメッシュ分割したモデル（混合体モデル）として構成されている。

混合法用腕モデル５０を用いて解析をおこなう場合、身体内部情報計算部１３は、腕モデル３０ａのＦＦＤ変形結果から得られる骨、筋肉、皮膚の位置を混合法用腕モデル５０に与える。この際、モーフィング等の手法を用いて、ＦＦＤから得られる骨、筋肉、皮膚の形状データを混合法用腕モデル５０に与える。身体内部情報計算部１３は、混合法用腕モデル５０に含まれる骨５１、筋肉５４、皮膚５３に対して、各時間ステップで通常の有限要素解析をおこなう。脂肪５２、神経５６、血管５５の混合体５７に対しては、以下で説明する三相理論による混合法を用いて解析する。三相理論による混合法としては、例えば、非特許文献２の技術を採用することができる。

三相理論による混合法では、固体相・流体相・イオン相が巨視的に同じ空間位置を同時に占有すると仮定し、各相間の力学的・電気的な相互作用を考慮した平衡方程式を用いる。ここでは、脂肪５２と神経５６と血管５５との混合体５７において、脂肪５２、血管５５の血管壁、神経５６の神経軸索を固体相、血管５５内の血液、神経５６内の神経伝達物質の流れを流体相、神経伝達物質が細胞へと伝わるときの過程をイオン相としてモデル化し、メッシュの各節点における平衡方程式から未知変数の解を求めることで解析をおこなう。

支配方程式は混合体・流体・イオンの３つの平衡方程式と３つの付帯条件式からなり、以下の式（４）〜（９）となる。式（４）は混合体の平衡方程式であり、式（５）は固体の圧縮条件式であり、式（６）は流体の平衡方程式であり、式（７）は混合体の圧縮条件式であり、式（８）はイオンの平衡方程式であり、式（９）は電気的中性条件である。ここではメッシュの各節点における未知変数である、固体相の変位（３自由度）をｕ、流体の相対速度（３自由度）をＱ^w、固体相の圧力（１自由度）をλ^S、混合体の圧力（１自由度）をλ^m、４種類のイオンと５種の代謝物質の濃度（９自由度）をｃα、電位（１自由度）をΨとしている。

ここで、χは変形前の固体相を基準とする参照配置である。Ｆは変形勾配である。Ｓは混合体の第２Piola-Kirchhoff応力である。Ｊは体積変化率である。κ^sは固体相の体積弾性率である。ｃαはイオンαのモル濃度である。ｃ^totはイオンモル濃度の合計である。ｃ^Fは単位流体体積あたりの固定電化である。ｚはイオンの価数である。Ψは静電ポテンシャルである。φ^wは混合体における流体の割合である。Φ^wは参照配置における体積に対する流体の割合（Φ^w＝Ｊφ^w）である。Ｄαはイオンαの拡散係数である。Ｋは透水係数（流体と固体の摩擦係数の逆数）である。φは浸透係数である。Ｒは気体定数である。Ｔは絶対温度である。Ｆｃはファラデー定数である。

これらの式（４）〜（９）から以下の式（１０）に示す支配方程式が導かれる。ただし、すべての変数について空間的に離散化し解析領域のすべての節点についてまとめた未知ベクトル｛Ｘ｝を｛Ｘ｝＝｛｛ｕ｝^T，｛Ｑ^w｝^T，｛λ^S｝^T，｛λ^m｝^T，｛ｃα｝^T，｛Ψ｝^T ｝^Tとする。

ここで、Ｑは内力ベクトルであり、Ｆは外力ベクトルである。この支配方程式を解くことで未知ベクトルＸを算出できる。しかし、力学・電位・イオンの反応はそれぞれ異なる時間スケールを持ち更新すべき頻度が異なるため、力学の機械系４変数（変位、流速、圧力のLagrange未定乗数２つ）と、電位と、イオンと、を分離すると、それぞれ以下の式（１１）ように表すことができる。

これらのベクトル｛Ｘ_m｝，｛Ｘ_p｝，｛Ｘ_c｝に対して、それぞれ異なる時間ステップΔｔｍ，Δｔｐ，Δｔｃを適用して更新する。イオンは現象の時間スケールが最も短く、電位はその次に小さい時間ステップであり、力学は最も長い時間スケールを持つ。そのため、電位やイオンについては各時間ステップにおいて常に解を更新する必要はない。このようにして、各時間ステップで有限要素解析をおこなうと、流体相でモデル化した血流速度や神経伝達物質の速度が得られ、そこから、血管血流量や神経伝達物質の移動量を得ることができる。

なお、神経軸索に跳躍伝導が生じる場合、神経伝達物質の流れる速度は血流などに比べて大きいため、混合法を用いるメリットが小さくなる。そのため、この場合には固体相と流体相のみを利用する二相理論による混合法を採用してもよい。ただし、神経伝達が神経軸索からシナプスを介して化学的伝達により次の細胞に移動する過程はイオンでモデル化できるため、この部分をモデル化する場合には、イオン相を含めた三相理論による混合法を用いる法が好ましい。これによりイオンの変化から電位場も算出できる。

以上説明した、第１実施形態の身体内部情報推定装置１０によれば、身体運動追従変形計算部１２が、剛体リンクモデル２０の身体運動に追従して身体内部組織の形状変化を算出し、身体内部情報計算部１３が、身体内部組織の形状変化からマクロ計算またはミクロ計算によって身体内部情報を推定するため、身体運動情報から身体内部情報を推定することができる。また、本実施形態の身体内部情報推定装置１０によれば、身体内部情報計算部１３は、数理モデルを用いて、身体内部情報のうちの血圧Ｐと筋力Ｆ^Mを算出するマクロ計算部１３１と、数値解析によって、身体内部情報のうちの血管血流量と神経電流量を算出するミクロ計算部１３２と、を含んでいるため、マクロ計算部１３１によって低計算コストで身体内部情報を算出できるとともに、ミクロ計算部１３２によって精確な身体内部情報を算出できる。

このような、本実施形態の身体内部情報推定装置１０によれば、身体の静止中または動作中において、既存の身体計測装置では計測が困難な身体内部情報を、コンピュータ上の人体モデルおよび人体モデルに身体機能を与える計算プログラムによって得ることができる。これにより、例えば、ゴルフスイング中において、関節における屈曲・伸展運動のともなう上腕の筋肉、血管、神経などの身体内部組織の形状変形と、その際の各組織の応力・ひずみ状態、血管・神経の伝達量、筋力等を算出できる。また、心拍運動や呼吸運動の時系列変化において、筋肉、血管、神経等の身体内部組織の形状変化と、その際の各組織の応力・ひずみ、血管・神経の伝達量、大動脈の血圧等を算出できる。

また、本実施形態の身体内部情報推定装置１０によれば、運動による筋負荷、関節負荷だけでなく、心拍上昇、血圧、呼吸変化等、身体運動中の身体全体にかかる負荷を可視化することができるので、スポーツ競技選手の育成強化や呼吸器系、循環器系を含む様々な疾患を抱える患者のリハビリ支援や、高齢者の歩行支援等に利用することができる。また、既存の身体計測装置では計測が困難な身体内部情報を低計算コストで算出できるので、繰り返し学習に活用しやすく、身体内部情報を取り入れたニューラルネットワークモデルを容易に構築することができる。

＜第２実施形態＞
図１２は、第２実施形態の身体内部情報計算部１３Ａの構成を例示した説明図である。第２実施形態の身体内部情報計算部１３Ａは、第１実施形態の身体内部情報計算部１３と比較すると、さらに、統合部１３３を備えている点が異なる。統合部１３３は、マクロ計算部１３１からの身体内部情報と、ミクロ計算部１３２からの身体内部情報を相互に補間しながら利用する。例えば、統合部１３３は、マクロ計算部１３１からの身体内部情報をミクロ計算部１３２に与え、ミクロ計算の境界条件に利用させることができる。一方、ミクロ計算部１３２からの身体内部情報をマクロ計算部１３１に与え、マクロ計算の式の物性パラメータに利用させることができる。

＜筋力の算出例＞
第１実施形態では、マクロ計算部１３１が、式（１）により筋力を算出する際に、筋断面積ＰＣＳＡ^Mは一定値と仮定した。しかし、第２実施形態では、マクロ計算部１３１が、第１実施形態の式（１）により筋力を算出する際に、統合部１３３は、ミクロ計算部１３２が算出した筋断面積Ａをマクロ計算部１３１に与える。すなわち、マクロ計算部１３１は、ミクロ計算部１３２が算出した筋断面積Ａを式（１）の筋断面積ＰＣＳＡ^Mとして使用する。ミクロ計算部１３２は、腕の筋断面積Ａを算出する場合には、混合法用腕モデル５０（図１１（ｂ））を用いて、各時間ステップにおいてＦＦＤ計算で得られる骨、筋、皮膚等の位置情報を、混合法用腕モデル５０の境界条件として与える。混合法用腕モデル５０の筋や皮膚を超弾性体などの材料でモデル化すれば、血液など流体相の影響を考慮して筋の変形を算出することができるので、より精度よく筋断面積Ａを求めることができる。また、混合法の有限要素解析で腕の変形状態を解析する場合には、ＦＦＤで腕を変形させた場合と異なり、変形時に筋が骨に貫入する事態が生じにくく、この意味でもより精度よく筋断面積Ａを算出できるといえる。

図１３は、羽状筋６４を例示した説明図である。式（１）の筋断面積ＰＣＳＡ^Mは、図１３に示すように、筋線維の方向の断面積であるため、羽状角θを考慮した断面積であるＡｃｏｓθがより好ましい。本実施形態によれば、筋の３次元的な詳細な変形状態から解剖学図の筋走行を参照して羽状角を算出できるため、より高精度に筋力を算出することができる。

＜力・応力・ひずみや血流量の算出例＞
ミクロ計算部１３２は、上述した混合法用腕モデル５０（図１１（ｂ））を用いた解析により、筋や骨にかかる力・応力・ひずみだけでなく、血管や神経にかかる力・応力・ひずみや血流量も計算することができる。この際、統合部１３３は、マクロ計算部１３１によるマクロ計算の結果をミクロ計算部１３２に与えてもよい。これにより、算出される値の精度の向上を図ることができる。統合部１３３は、マクロ計算部１３１によるマクロ計算から得られた骨や筋の力、血流量などの情報を、混合法用腕モデル５０の境界面における骨と筋の力や血流量（境界条件）としてミクロ計算部１３２に入力してもよい。また、骨と筋の形状の時系列変化をミクロ計算部１３２に与えれば、混合法有限要素解析の計算は、脂肪層と血管と神経の部分だけとなりミクロ計算部１３２の計算量を減らすことができる。

＜血圧の算出例＞
マクロ計算部１３１が、第１実施形態の式（２）により血圧を算出する際に、統合部１３３は、ミクロ計算部１３２が算出した大動脈のコンプライアンスＣｗと抹消血管抵抗Ｒｐをマクロ計算部１３１に与えてもよい。これにより、マクロ計算部１３１は、ミクロ計算部１３２が算出した大動脈のコンプライアンスＣｗを式（２）のコンプライアンスθとして使用することができる。第１実施形態では、式（２）において、大動脈のコンプライアンスθと抹消血管抵抗Ｒは定数としている。しかし、大動脈のコンプライアンスＣｗは、以下の式（１２）に示すように、最高血圧時の大動脈の断面積の最大値Ａ_max、動脈の効果長ｌ、および、心臓の収縮、拡張におけるそれぞれの近心動脈内血圧ｐの影響を受ける。

上記式（１２）の近心動脈内血圧ｐ、および、式（２）の末梢血管抵抗Ｒは、以下で説明するように、末梢での神経伝達物質ノルエピネフリン（ノルアドレナリン）の濃度Ｃｖ_NEの影響を受ける。具体的には、近心動脈内血圧ｐ、および、末梢血管抵抗Ｒは、下記の式（１３）、および、式（１４）に示すように、ともに式（１５）で定義される拡張期血圧を決定づける因子τｖ（ｔ）の関数になっている。式（１５）および式（１６）に示すように、この因子τｖ（ｔ）は、末梢での神経伝達物質ノルエピネフリンの濃度Ｃｖ_NEの影響を受ける。ここから、近心動脈内血圧ｐ、および、末梢血管抵抗Ｒは、末梢での神経伝達物質ノルエピネフリンの濃度Ｃｖ_NEの影響を受けることがわかる。

ここで、τｖ_NEは定数、ｋ^s _CvNEはは交感神経の活動度ｖ_sの重み、θ_vNEは時間遅れである。

以上から、大動脈の断面積Ａ、大動脈の効果長さｌ、神経伝達物質ノルエピネフリンの濃度Ｃｖ_NEがわかれば、末梢血管抵抗Ｒ、大動脈のコンプライアンスθが得られ、血圧をより精度よく算出できることがわかる。最高血圧時の大動脈の断面積の最大値Ａ_max、および、動脈の効果長ｌは、大動脈の伸展性や容積変化の特徴を示している。大動脈の伸展性が低下すると血圧が上昇し、大動脈が硬いと容積変化に対応するため血圧が上昇する。大動脈のコンプライアンスθは、単位圧力の増加に対する血管容積の増加量を示す。末梢血管抵抗Ｒは、血管の硬さと関係する。つまり、上記によって末梢血管抵抗Ｒ、および、大動脈のコンプライアンスθを得ることで、血管の硬さを考慮することができ、血圧の予測精度の向上を図ることができる。

ＦＦＤと混合法による３次元的な大動脈の変形計算によって、血管の収縮を考慮した、動脈の断面積Ａや動脈の効果長さlを算出できる。また、混合法によって神経伝達物質ノルエピネフリン濃度Ｃｖ_NEを算出できる。よって、ＦＦＤと混合法から、末梢血管抵抗Ｒと大動脈のコンプライアンスθを算出でき、精度のより高い動脈血圧Ｐを算出することができる。ノルエピネフリン濃度Ｃｖ_NEは、例えば、混合法用腕モデル５０において、神経伝達をシナプス伝達まで考慮してイオン相を扱うと算出することができる。なお、ノルエピネフリン濃度Ｃｖ_NEは、交感神経の活動度や時間遅れの影響を受ける。混合法を用いてノルエピネフリンのイオン濃度を算出し、ＦＦＤによる神経の変形から時間遅れを算出することで、交感神経の活動度の算出過程で仮定された定数を特定することができ、算出精度を向上させることができる。

以上説明した、第２実施形態の身体内部情報推定装置によれば、身体内部情報計算部１３Ａの統合部１３３は、マクロ計算部１３１からの身体内部情報と、ミクロ計算部１３２からの身体内部情報を相互に補間しながら利用するため、身体内部情報の精度を向上させつつ計算コストの低減を図ることができる。例えば、統合部１３３は、マクロ計算部１３１によって数理モデルを用いて算出された身体内部情報をミクロ計算部１３２に与え、ミクロ計算の境界条件に利用させることができる。これにより、数値解析によって算出される身体内部情報の計算コストを低減させることができる。一方、ミクロ計算部１３２によって数値解析によって算出された身体内部情報をマクロ計算部１３１に与え、マクロ計算の式の物性パラメータに利用させることができる。これにより、数理モデルを用いて算出される身体内部情報の精度を向上させることができる。

＜本実施形態の変形例＞
本発明は上記の実施形態に限られるものではなく、その要旨を逸脱しない範囲において種々の態様において実施することが可能であり、例えば次のような変形も可能である。

［変形例１］
第１、第２実施形態で示したマクロ計算部１３１によるミクロ計算例や、ミクロ計算部１３２によるマクロ計算例は、その一例であり、マクロ計算部１３１やミクロ計算部１３２によって推定される身体内部情報は、本実施形態で例示した内容に限定されない。例えば、非特許文献１に記載されている数理モデルを用いて身体内部情報を算出してもよいし、非特許文献２に記載されている数値解析法を採用してもよい。また、第２実施形態で示した統合部１３３は、マクロ計算部１３１によって算出された筋力や血圧と、ミクロ計算部１３２によって算出された筋肉の応力、ひずみ、血管の血流量、流速、圧力、神経の電流量、イオン伝達速度から、関節角度、筋伸び、血圧、心拍数等の身体内部情報を算出してもよい。

［変形例２］
第１実施形態では、混合法用モデルとして腕モデル５０を例示した。しかし、混合法用モデルは、身体の任意の部位について作成可能である。すなわち、身体内部情報計算部１３は、身体の任意の部位についての身体内部情報を推定することができる。また、本実施形態では人体モデルを用いて人間の身体内部情報を推定しているが、身体内部情報推定装置１０は、人以外の任意の生き物のモデルを用いて、その生き物の身体内部情報を推定することができる。

［変形例３］
第１実施形態の身体内部情報推定装置１０は、記憶装置に、剛体リンクモデル２０、人体モデル３０、筋骨格モデル４０、混合法用腕モデル５０が記憶されているものとして説明した。しかし、これらのモデルは例示であり、これらの一部のモデルを記憶していなくてもよいし、これら以外のモデルを記憶していてもよい。また、身体内部情報推定装置１０は、身体内部情報を算出する際に、必要な身体モデルをその都度作成する機能を有していてもよい。

［変形例４］
第１実施形態の身体内部情報計算部１３は、マクロ計算部１３１とミクロ計算部１３２を備えているものとした。しかし、身体内部情報計算部１３は、マクロ計算部１３１とミクロ計算部１３２の一方のみを備えていてもよい。

［変形例５］
第１実施形態の身体内部情報推定装置１０は、マクロ計算部１３１によるマクロ計算と、ミクロ計算部１３２によるミクロ計算を各タイムステップで同時におこなうこともできるし、マクロ計算のみ、ミクロ計算のみをおこなうこともできる。また、マクロ計算、ミクロ計算の同時計算について繰り返し計算をおこなうことや、身体の複数部位（上肢、下肢、内臓など）の結果を利用しあうこともできる。これにより、マクロ計算、ミクロ計算それぞれの計算精度の向上を図ることができる。

以上、実施形態、変形例に基づき本態様について説明してきたが、上記した態様の実施の形態は、本態様の理解を容易にするためのものであり、本態様を限定するものではない。本態様は、その趣旨並びに特許請求の範囲を逸脱することなく、変更、改良され得ると共に、本態様にはその等価物が含まれる。また、その技術的特徴が本明細書中に必須なものとして説明されていなければ、適宜、削除することができる。

１０…身体内部情報推定装置
１１…身体運動情報入力部
１２…身体運動追従変形計算部
１３…身体内部情報計算部
１３１…マクロ計算部
１３２…ミクロ計算部
１３３…統合部
１４…身体内部情報出力部
２０…剛体リンクモデル
２１…リンク
２２…ジョイント
３０…人体モデル
４０…筋骨格モデル
５０…混合法用腕モデル
３１、５１…骨
４１…筋肉モデル
３２…関節
５２…脂肪
３３、５３…皮膚
３４、５４…筋肉
５４１…上腕筋
５４２…上腕二頭筋
５４３…上腕三頭筋
３５、５５…血管
３６、５６…神経
３７…心臓
５７…混合体

Claims (10)

1. 骨格を表す剛体リンクモデルのデータと、
   前記剛体リンクモデルに関連づけられた、骨、筋肉、血管、および、心臓を含む身体内部組織のモデルによって構成される身体モデルのデータと、を用いて、血圧や筋力を含む身体内部情報を推定する身体内部情報推定方法であって、
   前記剛体リンクモデルに運動を表現させるための身体運動情報を入力する入力工程と、
   前記身体運動情報に基づく前記剛体リンクモデルの運動に追従して前記身体内部組織の形状変化を算出する変形計算工程と、
   数理モデルと数値解析の少なくとも一方を用いて、前記身体内部組織の形状変化から前記身体内部情報を推定する身体内部情報推定工程と、
   推定された前記身体内部情報を出力する出力工程と、を備え、
   前記変形計算工程では、前記剛体リンクモデルの運動に追従して前記身体モデルをＦＦＤ（Free Form Deformation）によって変形させることにより前記身体内部組織の形状変化を算出する、
   身体内部情報算出方法。
2. 請求項１に記載の身体内部情報推定方法において、
   前記身体内部情報推定工程には、
   数理モデルを用いて、前記身体内部情報のうちの血圧と筋力を算出するマクロ計算工程と、
   数値解析によって、前記身体内部情報のうちの血管血流量と神経電流量を算出するミクロ計算工程と、の少なくとも一方の工程が含まれている、
   身体内部情報推定方法。
3. 請求項２に記載の身体内部情報推定方法において、
   前記身体内部情報推定工程には、前記ミクロ計算工程において数値解析によって算出された前記身体内部情報を、前記マクロ計算工程において数理モデルを用いて前記身体内部情報を算出する際に利用する統合工程が含まれている、
   身体内部情報推定方法。
4. 請求項２に記載の身体内部情報推定方法において、
   前記身体内部情報推定工程には、前記マクロ計算工程において数値モデルを用いて算出された前記身体内部情報を、前記ミクロ計算工程において数値解析によって前記身体内部情報を算出する際の境界条件として利用する統合工程が含まれている、
   身体内部情報推定方法。
5. 請求項１から請求項４までのいずれか一項に記載の身体内部情報推定方法において、
   前記身体内部情報推定工程には、前記変形計算工程で計算された特定の筋肉の形状変化から、前記筋肉の長さの変化と、前記筋肉の収縮速度を算出し、算出した前記筋肉の長さの変化と、前記筋肉の収縮速度から数理モデルを用いて前記特定の筋肉の筋力を算出する工程が含まれている、
   身体内部情報推定方法。
6. 請求項１から請求項５までのいずれか一項に記載の身体内部情報推定方法において、
   前記身体内部情報推定工程には、前記変形計算工程で計算された心臓の形状変化から、左心室の容積変化を算出し、算出した前記左心室の容積変化から数理モデルを用いて血圧を算出する工程が含まれている、
   身体内部情報推定方法。
7. 請求項１から請求項６までのいずれか一項に記載の身体内部情報推定方法において、
   前記身体内部情報推定工程には、前記身体内部組織のうち、血管と脂肪とが一つの混合体として構成された混合体モデルを用いて、有限要素解析により血管血流量を推定する工程が含まれている、
   身体内部情報推定方法。
8. 骨格を表す剛体リンクモデルのデータと、
   前記剛体リンクモデルに関連づけられた、骨、筋肉、血管、および、心臓を含む身体内部組織のモデルによって構成される身体モデルのデータと、を用いて、コンピュータに血圧や筋力を含む身体内部情報を推定させるためのコンピュータプログラムであって、
   前記剛体リンクモデルに運動を表現させるための身体運動情報を入力する入力機能と、
   前記身体運動情報に基づく前記剛体リンクモデルの運動に追従して前記身体内部組織の形状変化を算出する変形計算機能と、
   数理モデルと数値解析の少なくとも一方を用いて、前記身体内部組織の形状変化から前記身体内部情報を推定する身体内部情報推定機能と、
   推定された前記身体内部情報を出力する出力機能と、をコンピュータに実行させるものであり、
   前記変形計算機能は、前記剛体リンクモデルの運動に追従して前記身体モデルをＦＦＤ（Free Form Deformation）によって変形させることにより前記身体内部組織の形状変化を算出する、
   コンピュータプログラム。
9. 記憶媒体であって、
   請求項８に記載のコンピュータプログラムを記憶する、
   記憶媒体。
10. 血圧や筋力を含む身体内部情報を推定する身体内部情報推定装置であって、
    骨格を表す剛体リンクモデルのデータと、前記剛体リンクモデルに関連づけられた、骨、筋肉、血管、および、心臓を含む身体内部組織のモデルによって構成される身体モデルのデータと、を記憶する記憶部と、
    前記剛体リンクモデルに運動を表現させるための身体運動情報を入力する身体運動情報入力部と、
    前記身体運動情報に基づく前記剛体リンクモデルの運動に追従して前記身体内部組織の形状変化を算出する変形計算部と、
    数理モデルと数値解析の少なくとも一方を用いて、前記身体内部組織の形状変化から前記身体内部情報を推定する身体内部情報推定部と、
    推定された前記身体内部情報を出力する出力部と、を備え、
    前記変形計算部は、前記剛体リンクモデルの運動に追従して前記身体モデルをＦＦＤ（Free Form Deformation）によって変形させることにより前記身体内部組織の形状変化を算出する、
    身体内部情報推定装置。