3712

チュートリアル

汎用マルチフィジックス有限要素解析ソフトウェア COMSOL Multiphysics を取り上げて新しい計算科学工学の仕事の流れを解説していきます。今回はチュートリアルの第4回目として、生物系を中心に解説していただきます。実験的な方法にも触れています。なお、前号よりチュートリアル記事は1ページ目のみを本誌掲載し、続きは日本計算工学会 HP 上で公開していますので、そちらも併せてご参照ください。

COMSOL Multiphysics による計算科学工学– 生物系（4）

橋口 真宜
三井 敏之
1

はじめに

間履歴を COMSOL Multiphysics のグローバル ODE イン
ターフェース（ODE: Ordinary Differential Equations）で初期値 x0（5、15、25、28）について計算しました。
どの初期値からも20という数値に漸近しています。
これは微分方程式が時間の発展形としてわかっていれば時間方向には図1に示した通り、xt-r*x*(1-x/k) と記述
すればすぐに計算できます。ただし、xt は x の時間 t に関する微分を表します。しかし、初期値 x0をもっと多くの値について検討したい場合にはこの方法では時間がかかってしまい、見通しも悪い方法です。

化学反応において可逆変化が化学種 A と B の時間に関する連立常微分方程式で記述されることはすでに見た 通 り で す。 平 衡 点 解 析 が Back-of-the-envelope
calculation として有効であることも見ました。この反応
を捕食者、被捕食者の関係としてみると、面白いことに生物や生態系の諸現象が化学反応と同じような考え方で記述できることになります。化学反応のところで空間依存モデルに触れました。生物や生態系でも拡散
項 に よ る 空 間 依 存 項 が 重 要 な 役 割 を 果 た し ま す。
チューリングの提唱した形態形成の仮説が実験でも検証されてきています。今回は生物や生態系について見ていきます。
1.1 ロジスティックモデルと平衡点解析
まず、マルサスが人口論で提唱した生物個体数が指数関数的に、あるいは幾何級数的に増殖するモデルに関する修正モデルであるフェアフルストのロジス
ティックモデルは次の形をしています [1]。r は成長率、
k は環境収容力です。r=0.5、k=20としたときの x(t) の時
筆者紹介
はしぐち まさのり
計測エンジニアリングシステム（株）第1技術部部

図1長。COMSOL Multiphysics の 端 緒 を 切 り 開 い た ス

ロジスティックモデルの時間解と平衡点解析

ウェーデンの数学者 Germund Dahlquist 先生とともに。COMSOL 社の受付ロビーにある先生の記念碑を出張の折に撮影し、ここに一緒に掲載しました。

続きは Web で
日本計算工学会誌「計算工学（Vol.23, No.1）」HP:
http://www.jsces.org/activity/journal/

第1章
みつい としゆき
1999年ミネソタ大学大学院物理学 Ph.D 修了、カリ
フォルニア大学、ハーバード大学博士研究員を経て現在、青山学院大学理工学部物理数理学科教員、
専門は実験物理学、特にナノデバイス、生物物理。
現象の解明に COMSOL Multiphysics を用いる。

第2章〜4章注．COMSOL, COMSOL Multiphysics, COMSOL Desktopおよび Livelinkは
COMSOL ABの商標または登録商標です。その他の製品名、サービ
ス名、組織名は、各組織の商標または登録商標です。なお、本文中に™、®マークは明記しておりません。本稿の説明は著者独自の考えに基づいており他を代表するものではありません。

（22）

計算工学

チュートリアルCOMSOL Multiphysics による計算科学工学 – 生物系（4）そこで、平衡点解析を考えてみます。これは dx/dt=0
チューリングは拡散方程式でも不安定化する場合があとなる点を検討する方法であり、x=0と x=k=20が平衡ることを示しました。その考え方は生物系の実験でも検点（定常解）であることがすぐにわかります。さらに、証されてきています。数値計算も COMSOL Multiphysics
v=dx/dt を ベ ク ト ル 場 と 考 え る と、 図1の よ う に、を 利 用 し て 数 多 く 行 わ れ て き て い ま す。COMSOL
r*x*(1-x/k) の符号によって +1、-1をとる折れ線グラフ
Multiphysics を活用した取り組みは DNA といった生物物
チュートリアルCOMSOL Multiphysics による計算科学工学 –生物系 (4)と一緒に描画すると、その符号によってベクトル場の理の領域でも広範囲に行われています。
向きが右を向くか左を向くかが決まっていることがわ
は周囲のベクトル場はその点を向いているので安定点の挙動解明への取り組み
かります。x=0に初期値をとったとすると、ここからは
2 ２．DNADNA の挙動解明への取り組みになり、沈点（シンク）です。つまり、時間積分を実
本章では、私、三井の行っている三つの『実験』と、常に x を増やす右方向にベクトル場が作用するので、本章では、私、三井の行っている三つの『実験』と、行しなくても単にグラフを描くだけで解が初期位置か
そのシミュレーションについて説明します。①ナノポ
湧点（ソース）と呼ばれる不安定点であることがわかりそのシミュレーションについて説明します。①ナノポ
らどの位置へ移動していくかがわかることになります。 アによる DNA 検出、②羽の原基形成と Turing の反応拡ます。一方で
x=k=20は周囲のベクトル場はその点を向アによる DNA 検出、②羽の原基形成と Turing の反応拡
これは大変有効な Back-of-the-envelope calculation
散方程式、③心筋細胞の拍動とそのシグナル伝搬、のいているので安定点になり、沈点
（シンク）
です。つま
散方程式、③心筋細胞の拍動とそのシグナル伝搬、のです。
三つです。実験を無駄なく効率的に行い、更に、現象り、時間積分を実行しなくても単にグラフを描くだけ
次のモデルを見てみましょう。遺伝子選択モデルにさ 三つです。実験を無駄なく効率的に行い、更に、現象を物理的に理解するために、計測が必須の時代になりで解が初期位置からどの位置へ移動していくかがわからに拡散モデルを合成したフィッシャーモデルがあり、 を物理的に理解するために、計測が必須の時代になりました。計測の教育は受けておらず、計測の詳細や、それは進行波解(travellingwave
ることに な り ま す。 こ れ は 大 変
有 効solution)を持ちまな Back-of-theました。計測の教育は受けておらず、計測の詳細や、それらの表現が適切ではないところもあるかと思います[1]。詳細は参考図書を参照していただくとして、結 それらの表現が適切ではないところもあるかと思いま
envelope
calculation です。
すが、ご了承をお願いいたします。それでは、まずは果の式を図２中に示します。
「 微分方程式は解かずに」、 すが、ご了承をお願いいたします。それでは、まずは次のモデルを見てみましょう。遺伝子選択モデルにナノポアについて説明します。
dP/dz、dQ/dz の右辺の関数をベクトル場の成分として ナノポアについて説明します。
さらに拡散モデルを合成したフィッシャーモデルがあ
その包絡線群を、COMSOL Multiphysics の流線描画機
t = 0.07
り、それは進行波解（travelling wave solution）を持ちま能を利用して(P,Q)平面に作図した結果も図２に示し
す [1]。詳細は参考図書を参照していただくとして、結ました。その向きは矢印を重ね書きしました。D=8、果の式を図2中に示します。
「微分方程式は解かずに」、
r=0.8、k=100、m=1 とした描画であり、c はコルモゴロ
チュートリアルCOMSOL Multiphysics による計算科学工学 –生物系 (4)dP/dz、dQ/dzの右辺の関数をベクトル場の成分としてフの速度にとっています。これからわかることは平衡
その包絡線群を、COMSOLMultiphysicsの流線描画機点(0,0)は安定結節点(Stablenode)には流線が周囲かは周囲のベクトル場はその点を向いているので安定点
２．DNA の挙動解明への取り組み
能を利用して(P,Q)
平面に作図した結果も図2に示しまになり、沈点（シンク）です。つまり、時間積分を実
本章では、私、三井の行っている三つの『実験』と、し た。
そ の 向 き は 矢 印 を 重 ね 書 き し ま し た。D=8、行しなくても単にグラフを描くだけで解が初期位置か
そのシミュレーションについて説明します。①ナノポ
r=0.8、k=100、m=1とした描画であり、cはコルモゴロ アによる DNA 検出、②羽の原基形成と Turing の反応拡らどの位置へ移動していくかがわかることになります。
これは大変有効な Back-of-the-envelope calculation
フの速度にとっています。これからわかることは平衡
散方程式、③心筋細胞の拍動とそのシグナル伝搬、のです。
点 (0,0)の安定結節点（Stable node）には流線が周囲から 三つです。実験を無駄なく効率的に行い、更に、現象次のモデルを見てみましょう。遺伝子選択モデルにさ流入し、他方の平衡点 (1,0) では鞍点（Saddle）になって を物理的に理解するために、計測が必須の時代になりらに拡散モデルを合成したフィッシャーモデルがあり、いることです。その平衡点同士を結ぶ軌道が流線の集ました。計測の教育は受けておらず、計測の詳細や、それは進行波解(travelling wave solution)を持ちま
それらの表現が適切ではないところもあるかと思いま
積として観察されます。これがフロント進行波解に相す[1]。詳細は参考図書を参照していただくとして、結図3 図３.
ナ ノ ポナノポアの実験模式図(左上)。ナノポア
ア の 実 験 模 式 図（左 上）。 ナ ノ ポ ア TEM 像（右[1]
すが、ご了承をお願いいたします。それでは、まずは上）。蛍光顕微鏡によるDNA 一分子観測。中心のポア
当するヘテロクリニック軌道です。これは理論解析
果の式を図２中に示します。
「 微分方程式は解かずに」、
TEM像(右上)。蛍光顕微鏡によるDNA 一分子に進入する様子（下図）。DNAの軌跡、中心にポア（右ナノポアについて説明します。
dP/dz、
dQ/dzの右辺の関数をベクトル場の成分としてと一致します。初期値を変化させて時間方向に計算す観測。中心のポアに進入する様子(下図)。DNA中）。DNA は電気泳動で進入する。
その包絡線群を、COMSOL Multiphysics の流線描画機の軌跡、中心にポア(右中)。DNA
は電気泳動で
t = 0.07
る方式ではこのような図を描くのに相当な計算量にな能を利用して(P,Q)平面に作図した結果も図２に示し進入する。
るでしょう。Back-of-the-envelope calculation として活用ました。その向きは矢印を重ね書きしました。D=8、
2.1
DNA
すべきです。数理物理学の考え方は有用です。
２. 実験：ナノポアによる
1 実験: ナノポアによる
DNA検出
検出
r=0.8、k=100、m=1 とした描画であり、c はコルモゴロ
半導体プロセスの発展により、マイクロサイズの加微分方程式を解かずに大域的な様子を知ることが可
半導体プロセスの発展により、マイクロサイズの加フの速度にとっています。これからわかることは平衡工が実験室レベルでも行われるようになり、材料も半能な場合があるということです。
点(0,0)は安定結節点(Stable node)には流線が周囲か工が実験室レベルでも行われるようになり、材料も半10 μm
導系の固体だけではなく、ガラスやポリジメチルシロ
導系の固体だけではなく、ガラスやポリジメチルシロ
図２ フィッシャーモデルの図式解法
キサン（PDMS）などのソフトマテリアルによるマイク
キサン（PDMS）などのソフトマテリアルによるマイク
ロサイズの加工が可能になりました。そして、流路やロサイズの加工が可能になりました。そして、流路やら流入し、他方の平衡点(1,0)では鞍点(Saddle)になっ容器の作製により、小さな化学工場としてのlab-on-a容器の作製により、小さな化学工場としてのていることです。その平衡点同士を結ぶ軌道が流線の
chip
や、 Nano-Micro-ﬂuidics によるセンサーデバイス“lab-on-a-chip”や、”NanoMicro-fluidics” に集積として観察されます。これがフロント進行波解に[2]1の開発が盛んになりました。また、単一細胞と同じよるセンサーデバイスの開発が盛んになりました。
相当するヘテロクリニック軌道です。これは理論解析
スケールの構造や流路により、単一細胞の仕分けや局
また、単一細胞と同じスケールの構造や流路により、[1]と一致します。初期値を変化させて時間方向に計算
ナノポアの実験模式 図(左上)。ナノポア
所的な刺激を可能となり、ミクロな細胞の研究に貢献
単一細胞の仕分けや局所的な刺激を可能となり、ミクする方式ではこのような図を描くのに相当な計算量に 図３.
TEM 像(右上)。蛍光顕微鏡による
DNA 2一分子[3]。現在では、医療診断から
Drug Delivery
ロな細胞の研究に貢献しています。現在では、医療なるでしょう。Back-of-the-envelope calculation と しています
観測。中心のポアに進入する様子(下図)。DNA
まで、分野を跨いだ領域で発展を遂げています。
診断から Drug Delivery まで、分野を跨いだ領域で発して活用すべきです。数理物理学の考え方は有用です。 の軌跡、中心にポア(右中)。DNA
は電気泳動で微分方程式を解かずに大域的な様子を知ることが可能 進入する。
マイクロ・ナノサイズの流路では、体積と表面積の展を遂げています。
な場合があるということです。
マイクロ・ナノサイズの流路では、体積と表面積のアスペクト比から、表面による効果が顕著に現れま２. 1 実験: ナノポアによる DNA 検出[4]チューリングは拡散方程式でも不安定化する場合があ す。その代表的な例が電気浸透流です
アスペクト比から、表面による効果が顕著に現れます。
。圧力差がな半導体プロセスの発展により、マイクロサイズの加ることを示しました。その考え方は生物系の実験でも くても流れが発生します。現在ではこの流れを簡易的
3
その代表的な例が電気浸透流です。圧力差がなくて工が実験室レベルでも行われるようになり、材料も半検 証 さ れ て き て い ま す 。 数 値 計 算 も COMSOL 10 μm[5]も流れが発生します。現在ではこの流れを簡易的なポなポンプとして応用するデバイスもあるくらいです。
導系の固体だけではなく、ガラスやポリジメチルシロ
図2 フィッシャーモデルの図式解法
Multiphysicsを利用して数多く行われてきています。
図２ フィッシャーモデルの図式解法
ンプとして応用するデバイスもあるくらいです 4 。
キサン（PDMS）などのソフトマテリアルによるマイク
COMSOL Multiphysics を活用した取り組みは DNA とい
私はナノポアを用いた DNA センサーの基礎研究をし
ロサイズの加工が可能になりました。そして、流路や計算工学
Vol.23, No.1 2018
った生物物理の領域でも広範囲に行われています。
ら流入し、他方の平衡点(1,0)では鞍点(Saddle)になっ （22-2）ていることです。その平衡点同士を結ぶ軌道が流線の集積として観察されます。これがフロント進行波解に計算工学相当するヘテロクリニック軌道です。これは理論解析[1]と一致します。初期値を変化させて時間方向に計算容器の作製により、小さな化学工場としての“lab-on-a-chip”や、”Nano- Micro-fluidics” に(??-2)
Vol.22,
1
よるセンサーデバイスの開発が盛んになりました。 No.4 2017

また、単一細胞と同じスケールの構造や流路により、

チュートリアルCOMSOL Multiphysics による計算科学工学 – 生物系（4）

私はナノポアを用いた DNA センサーの基礎研究をしています（図3）[6]。このナノポアデバイスの最終ゴール
は DNA シークエンスです。ナノポアはシリコン化合物
によるフリースタンディング薄膜（15~200nm）をつく
り、その膜にナノサイズの穴（ナノポア）をあけます（図3左上）。このナノポアに電気泳動させた DNA を通して、一分子レベルによる DNA の解析をします（図3下）[7]。ポアを金コートして、その電位を gate のように変化させて、DNA の挙動制御も試みています [8]。これまでのナノポア +DNA の研究結果から、シークエンスデバ
イスを完成させるためには、ナノスケールにおけるDNA とポアの壁（シリコン化合物や金）との相互作用について知る必要があります。理由はポア付近で、DNAが予期せぬ動きをするからです。DNA がポアに停滞したり、DNA 同士がポア付近で反発し合ったり、これらの現象は、ポア付近の電場による電気泳動だけではな
く、ポアや DNA を介した電気浸透流の影響も考えられます。実験的には、蛍光分子を DNA に直接つけて、光学顕微鏡を用いて DNA のダイナミクスを観測します。
しかし、空間的にも時間的にも解像度に限界があり、電場 と 流れ を分けて見積もれません。そこで、COMSOL Multiphysics を用いて、ナノポア付近の物理を調べる必要があるのです [8]。ナノポアの研究グループでは COMSOL がシミュレーションの代名詞で、近年ではマイクロポアを通過するバクテリア・ウイルス検出のシミュレーションにも用いられています [9]。
2.2 シミュレーション：DNA のポア通過とその物理学
ナノポア付近では、K+ や Cl- イオンの濃度に分布が生じで電場が形成されます。この電場に負の電荷を持つ
DNA が引かれ、ポアを通過します。ポア付近ではポア
表面の電荷や、金コートの場合はその電位により、電気浸透流の影響も予測できます。更に DNA 自体の電荷によっても流れが発生する可能性も有ります。そこ
で、定量的にポア付近の電位、流れを知る必要があります。Multiphysics な COMSOL の力を借りましょう。
一般的にマイクロ・ナノデバイスに共通する物理学として、イオン密度や電極由来の電場（Poisson Boltzmann
equation）、イオンの拡散や電場による移動等（NernstPlanck equation）、イオンの移動が引き起こす水溶液の流れ（Navier-Stokes equation）です。
この三つの式を iteration として定常的な解を求めます。その結果から、溶液内の電場と流れがわかり、電荷を持つ DNA などの分子や粒子のダイナミクスを見積もることができます。それぞれを簡単に説明します。
専門家の方は飛ばしてください。
① イ オ ン 密 度 や 電 極 電 位 に よ る 電 場（Poisson
Boltzmann equation）
ρe：電荷密度、φ：ポテンシャル、εr：比誘電率、
ε0：真空電率
イオン溶液中の電荷密度はイオンの濃度分布です。
極性に注意して、マイナスイオンの場合は計算工学

c：イオンの密度、e：電気素量、z：価数、
i：イオン種 (K+, Cl-)とします。シミュレーションでは価数に極性をつけます。
②

イオンの移動（Nernst-Planck equation）
電荷を持つイオンや分子等の溶液中の泳動を記述す
る電気泳動（electrophoresis）の式で、

μm：電気泳動移動度、D：イオンの拡散係数、
u：溶液の速度自由拡散にイオンの電気泳動の項を追加しただけです。厳密には移動度を含む項はイオンの極性により符号が変わります。①と②の式で、固体表面付近の水溶
液中のイオンが、固体の表面電荷（絶縁体）や、印加電位（導体）と溶液との電位差を遮蔽する電気二重層を見積もります。この二重層の厚さは、シンプルな形状を仮定すれば、解析的にも解けて、 デバイ長 として10-2M
程度のイオン濃度では3nm 程度とわかっています。
水溶液中のイオン種が K+ と Cl- のみの場合、それぞれの移動度の大きさが等しく、極性が逆のため、電場によるイオンの泳動で水溶液の流れは起きません。つまり、u は0です。しかし、遮蔽イオンの領域では偏りが出来て、流れ（電気浸透流）が発生します。この影響
は流体力学で記述します。
③

イ オ ン の 移 動 が 引 き 起 こ す 流 れ（Navier-Stokes
equation）

μ：粘性、ρ：密度、p：圧力
シリコン化合物では表面が負に帯電するので、遮蔽は + イオンです。このイオンによる電気浸透流が、この Navier-Stokes 式で見積もれます。pH7~8にて−の電荷となる DNA の電気泳動とは逆向きの流れとなります。ナノポアのグループは、この逆流を正確に見積もる必要があります。あるいは表面電荷を0にするコー
ティングを行うか。シミュレーションでは、①、②、③において、イオンの電荷密度、ポテンシャル、そして溶液の流れが共通の変数です。ちなみにポアを長いトンネル（channel）、デバイ長 << ポアの径と仮定すると、電気浸透流を、ポア方向のみの r による変形ベッセ
ル関数として見積もれます [10]。しかし、実際の実験条
件を考慮すると、自由度を削りたくないものです。
2.3 チュートリアル：DNA のポア通過とその物理学紹介した三つの式はそれぞれ、COMSOL では、AC/DC の静電場、化学種輸送の希釈種輸送、流体流れの層流を用います。本来、収束する定常的な解を求めた
く、①、②の式のみなら定常解が求まりますが、ナノ（22-3）

Vol.23, No.1 2018

気浸透流を、ポア方向のみの r による変形ベッセル関

定義して、共通の関数として用いると便利です。また、数として見積もれます 7 。しかし、実際の実験条件を∫Integration の関数もグローバル定義して、図４の考慮すると、自由度を削りたくないものです。

左中段にあるジオメトリの線 5 を測定用にとり、この線（ポア断面）を通るイオンの flux を積分して計測中

２．３ チュートリアル：DNA のポア通過とその物理
チュートリアル学に確認します。この値は実験で測定するイオン電流でCOMSOL Multiphysics による計算科学工学 – 生物系（4）す。図３の模式図の中の A と表した電流に相当します。
COMSOL で は そ れ ぞ れ の 方 程 式 ( ス タ デ ィ ) に ク ラ ス

space_charge の定義

図4

図４. 三つの式と COMSOL。N-P eq., N-S eq.を示し三つの式と COMSOL。N-P eq., N-S eq. を示しましました。space_chargeの関数を電荷密度として定義た。space_charge の関数を電荷密度として定義し依存として計算する、あるいは電位差を与える電極表
面やポア表面の電荷を、t=0では0として、sigmoid 関数
義の K イオン。流速まで存在する)として、細かく値がのように0から連続的に変化させていく方が、time step
貯めこまられており、簡単に値を拾えます。円柱座標
が無難に進みます。私はCOMSOL 上では step 関数があ系なので∫の際に2πr と電荷の正負はよく間違いまり、それに Smoothing をかけることができるので、t=0す。このイオン電流の値が計測時間の経過と共に一定で値が0となるように
step の場所をシフトして使っていの値に落ち着いたら、COMSOL スタディの定常解を待たます。
ずに、時間依存のまま終えてよい指標となります。
計測時間短縮に関しては、ポアの表面電荷や形状に今回、説明はしませんが、金コートの電位をシミュ
よりますが、静電場と希釈種輸送だけの定常解を求め
レーションに反映させるのには苦労しました。金表面て、その解を初期条件として層流を含む三つ式を時間の電位を実験パラメータと一致させるような表面電荷依存として計算する、あるいは電位差を与える電極表
分布を、常に見積もりながら、 表面電荷 として取り面やポア表面の電荷を、t = 0 では 0 として、sigmoid
扱わないと解までたどり着けませんでした。しかし、関数のように 0 から連続的に変化させていく方が、
そ ん な ト リ ッ ク も COMSOL は 可 能 に し ま す。 グtime step が無難に進みます。私は COMSOL 上では step
ローバル ODE で電位に対応する表面電荷を計算させチュートリアル関数があり、それに Smoothing をかけることができるて、その表面電荷分布を計算に取り込ませました。
(例: chds.fluxz_cK。z は z 成分で、cK はユーザー定COMSOL Mu

ので、t = 0 で値が 0 となるように step の場所をシフトして使っています。
①
絶縁体膜ポア
2.4 実験結果との考察：DNA のポア通過とその物理学
今回、説明はしませんが、金コートの電位をシミュします。このN-Seq.
eq.のにクーロン力としてます。この N-Sの FFにクーロン力として
space_
ここでは、絶縁体膜ポア、金コートのポア、そしレーションに反映させるのには苦労しました。金表面
charge∇V。
spacecharge∇V。
て、DNA がポアに停滞中のポアの三つについて説明しの電位を実験パラメータと一致させるような表面電荷
サイズの構造に、それを分割するメッシュで計算する 分布を、常に見積もりながら、“表面電荷”として取ります。
チュートリアル
紹介した三つの式はそれぞれ、COMSOL では、AC/DC
層 流（Navier-Stokes
equation）は 落 ち 着 つ き ま せ ん 。 扱わないと解までたどり着けませんでした。しかし、① 絶縁体膜ポアの静電場、化学種輸送の希釈種輸送、流体流れの層流
COMSOL は可能にします。グロー
そこで、時間依存（例えば t=10-3まで）として見切りを そんな“トリック”もの電気浸透流よる泳動（vosm）< 電気泳動（vel）と①DNA
絶縁体膜ポアを用います。本来、収束する定常的な解を求めたく、
バル
ODE
で電位に対応する表面電荷を計算させて、そ付けます。主な COMSOL のポイントとなる画面を切り①、②の式のみなら定常解が求まりますが、ナノサイ
出 し て 図4に ま と め ま し た。 先 ず は 左 上 の 三 つ の 式 の表面電荷分布を計算に取り込ませました。
ズの構造に、それを分割するメッシュで計算する層流（フィジクス）の従属変数を統一します。ジオメトリは(Navier-Stokes equation)は“落ち着つきません”。そ
２．４ 実験結果との考察：DNA のポア通過とその物理
円 こで、時間依存(例えば
柱 座 標 系 に て θ 依 存tな（2D
軸 対 称）と し ま す。
-3=し10まで)として見切りを学
メッシュは図のように最適化します。COMSOLは計算付けます。主な COMSOL のポイントとなる画面を切り出
ここでは、絶縁体膜ポア、金コートのポア、そして、
図５. 絶縁膜ポア付近の電場。DNA の挙動による電する式（スタディ）や構造のスケール、測定のための境して図４にまとめました。
先ずは左上の三つの式(フィ

さて COMSOL の

気浸透流よる泳予想される DNAは、電気浸透流

COMSOL Mul

縁体膜ポアの結した。
ii. Vのgat
さて COMSOL

から離れる
DNA
気浸透流よる泳動予想される DNA のは軸対称のシミ

は、電気浸透流がきを予測しまし

縁体膜ポアの結果
ポアの限られたした。 ii. V

gateに立ち上がり、から離れる DNA のの対称性のずれ
DNA がポアに停滞中のポアの三つについて説明します。
は軸対称のシミュ界線によりメッシュサイズを場所により最適化しまジクス)の従属変数を統一します。ジオメトリは円柱座
場と COMSOL による見積もりが一致。
る動き(流れ)がきを予測しました
標系にてθ依存なし(2Dす。
私 は 無 理 に マ ニ ュ軸対称)とします。メッシュは
ア ル で サ イ ズ を 変 え ず に、は3D でのシミュ
ポアの限られた断図のように最適化します。COMSOL
は計算する式(スタCOMSOL の推奨する要素サイズや分布に従い、同じ数 DNA の電気浸透流よる泳動( v osm )＜電気泳動( v el )となに立ち上がり、還ディ)や構造のスケール、測定のための境界線によりメ（factor）
をかけて変える程度にします。ただし、膜表面
り、ポア付近のわずかな電場の染みだしに
が引か
図5図５.
絶 縁絶縁膜ポア付近の電場。
膜 ポ ア 付 近 の 電 場。DNAの の挙動による電挙 動DNAによる電場とDNAッシュサイズを場所により最適化します。私は無理に③
DNA がポアにの対称性のずれ(ではメッシュサイズを nm 以下にして、Stern layer 的な れます(図５)。その速度から電場を見積もると、実験
COMSOLによる見積もりが一致。
場とCOMSOLによる見積もりが一致。
マニュアルでサイズを変えずに、COMSOL の推奨する要る動き(流れ)が観
ポテンシャルを得たいので、境界層メッシュ（四辺形）結果と、COMSOL の電位のプロファイルが一致します 5b 。
なり、ポア付近のわずかな電場の染みだしに DNA が引 は 3D でのシミュにします（図4）。私自身メッシュのとりかたにはよく悩
DNA の電気浸透流よる泳動( v osm )＜電気泳動( v el )とな
表面電荷の影響が限定的なので、球対称となれば、点かれます（図5）。その速度から電場を見積もると、実験[11]
みます。特筆すべき点は、space_charge
を図4のようり、ポア付近のわずかな電場の染みだしにDNA が引か
計算工学(??-4)
Vol.22, No.4 2017
電荷の電位とほぼ一致します。
結果と、COMSOL の電位のプロファイルが一致します③ DNA がポアににグローバル定義して、共通の関数として用いると便れます(図５)。その速度から電場を見積もると、実験[7]。表面電荷の影響が限定的なので、球対称となれば、
結果と、COMSOL の電位のプロファイルが一致します 5b 。
利です。また、∫ Integration の関数もグローバル定義し点電荷の電位とほぼ一致します。
②
金コートのポア
表面電荷の影響が限定的なので、球対称となれば、点て、図4の左中段にあるジオメトリの線5を測定用にと電荷の電位とほぼ一致します。
り、この線（ポア断面）を通るイオンの ﬂux を積分して 図３は金を挟んで膜をつくったポアの例ですが、この② 金コートのポア
計測中に確認します。この値は実験で測定するイオン ポアは、片面金の SiN-Au 膜です。実験結果は、V gate
図3は金を挟んで膜をつくったポアの例ですが、この②
金コートのポア
FET の gateV と考えています。実験
電流です。図3の模式図の中の A と表した電流に相当し はポア膜の電位でポアは、片面金のSiN-Au膜です。実験結果は、Vgate はポ
図３は金を挟んで膜をつくったポアの例ですが、この結果は、ます。COMSOL ではそれぞれの方程式（スタディ）にク
ア膜の電位でFETのgateVと考えています。実験結果は、
ポアは、片面金の SiN-Au 膜です。実験結果は、V gate
ラ ス（例：chds.ﬂuxz_cK。z は z 成 分 で、cK は ユ ー ザ ーi.
Vgate=V 溶液FET の gateV と考えています。実験
V
< V溶液
溶液
Vgate
gate<V
はポア膜の電位で27.94 s定義の K イオン。流速まで存在する）として、細かく値 i.
V gate = V 溶液
絶縁体膜ポアの結果と同じ
結果は、が貯めこまられており、簡単に値を拾えます。円柱座
絶縁体膜ポアの結果と同じ
ii.
Vgate<V 溶液
標系なので∫の際に2π r と電荷の正負はよく間違いま ii.
V gate < V 溶 液
V27.94 s2
gate < V 溶液i.
V gate =DNA
V 溶液ポアに入る（赤）中心付近のす。このイオン電流の値が計測時間の経過と共に一定
中心付近の
DNA
ポアに入る(赤)それ以外のDNAはポアから離れる（緑）
絶縁体膜ポアの結果と同じの値に落ち着いたら、COMSOL スタディの定常解を待
それ以外の
DNA
iii.
>V
ii. VVgate< 溶液
Vはポアから離れる(緑)gate溶液
図７. COMSOL
たずに、時間依存のまま終えてよい指標となります。
中心付近の
iii.
V gate> V 溶液 DNA ポアに入る(赤)6秒間ほどポアに入って、次の6秒間はポアから離れDNA を長くた
計測時間短縮に関しては、ポアの表面電荷や形状にそれ以外の DNA はポアから離れる(緑)る。ポアの上を一方向に流れる異方性のある動き
6 秒間ほどポアに入って、次の6 秒間はポアか
図７.
COMSOL 予
引き起こされるよりますが、静電場と希釈種輸送だけの定常解を求めiii.
V gate > V 溶液
DNA を長くたくら離れる。ポアの上を一方向に流れる異方性のて、その解を初期条件として層流を含む三つ式を時間6 秒間ほどポアに入って、次の 6 秒間はポアか引き起こされる実験では DNAある動きら離れる。ポアの上を一方向に流れる異方性の計算工学（22-4）

ある動き

Vol.23, No.1 2018

アに詰まってい実験では DNA がら DNA がポアにアに詰まっているもあれば、先にら DNA がポアに近うに観える場合
もあれば、先に停浸透流が実験的
うに観える場合も

それ以外の DNA はポアから離れる(緑)iii.

図７. COMSOL 予想と実験結果。ポアを大きく、

V gate > V 溶液

DNA を長くたくさんポアに置くと、DNA の逆流が6 秒間ほどポアに入って、次の 6 秒間はポアか引き起こされる。実験結果はそれを再現した。

ら離れる。ポアの上を一方向に流れる異方性のチュートリアルある動き

COMSOL Multiphysics による計算科学工学 – 生物系（4）
実験では DNA がポアに停滞することがあります、ポアに詰まっているかもしれません。そんな状況で後か

COMSOL では、ポアを大きく、DNA を 多く ポア
ら DNA がポアに近づくと、すり抜けてポアに入る場合に進入させると、DNADNA
由来の電気浸透流で、他のDNAもあれば、先に停滞しているに“はじかれる”よがポアに進入しない予想をしました。そして実
うに観える場合もあります。数年前にDNA （図7上）由来の電気

9
浸透流が実験的に観測されました。ちなみにその論
験で試しました。T4という購入できる最長（56.4μm）

文も COMSOL
そこで、私は“はじく”の DNAによる解析でした。
を、そしてポアの径を200nmにしました。実験
条件を、実験前に
COMSOL で予測しました。ポア内に結果ではポアに進入しようとした
T4DNA がポアにて一DNA を仮想的に置いて、電気浸透流の流れを観て、DNA
時停滞して、逆に手前に出てくる現象が観測されましが“はじかれる”ものでしょうか。
た（図7下）！ポアの径を大きくして初めて観測された
COMSOLでは、ポアを大きく、DNA を“多く”ポアに現象です。驚きました。
進入させると、DNA由来の電気浸透流で、他のDNA でのが DNA の
ナノポアの実験による研究ではCOMSOL

ポアに進入しない予想をしました（図７上）。そして実進入頻度や電流の見積もりは、default です。COMSOL
図6（青 vel）（緑 vosm）の和を
DNA の動図６.電気泳動
電気泳動(青v elと電気浸透流)と電気浸透流(緑
v osm
)の和を
験で試しました。T4 という購入できる最長(56.4 μm)
き（赤）として COMSOL の見積もり。
なくして研究をすすめられません。現在、DNA の異方
DNA の動き(赤)として COMSOL の見積もり。
の DNA を、そしてポアの径を 200 nm にしました。実
COMSOL Multiphysics による計算科学工学 –生物系 (4)さて COMSOL の結果を図6にまとめました。DNA の
さて COMSOL の結果を図６にまとめました。
DNA の電電気浸透流よる泳動（vosm）と電気泳動（vel）の速度合成
計算工学(??-5)
気浸透流よる泳動(
v osm )と電気泳動( v el )の速度合成をを予想される
DNA の動きとして、まず
i.Vgate=V 溶液 で予想される DNA の動きとして、まず i. V gate = V 溶液 では、電気浸透流が観測されなかったので、予想通り絶は、電気浸透流が観測されなかったので、予想通り絶縁体膜ポアの結果と同じになって実験結果と一致しま
縁体膜ポアの結果と同じになって実験結果と一致しました。ii.Vgate<V 溶液では、中心から外れるとポアから離した。 ii. V gate < V 溶 液 では、中心から外れるとポアれるから離れる
DNA の挙動が実験観測と一致しました。iii.
は軸対DNA の挙動が実験観測と一致しました。iii.
称のシミュレーションですが、大きな還流の動きを予
は軸対称のシミュレーションですが、大きな還流の動測しました。膜表面からポアに向かう浸透流が、ポアきを予測しました。膜表面からポアに向かう浸透流が、ポアの限られた断面積内に入りきれないvolumeが上の限られた断面積内に入りきれないvolumeが上に立ちに立ち上がり、還流を引き起こします。実際は、ポア
上がり、還流を引き起こします。実際は、ポアの対称
NA の挙動による電の対称性のずれ(楕円、角度)から、大きな異方性のあ性のずれ（楕円、角度）から、大きな異方性のある動きが一致。
る動き(流れ)が観測されたのだと思います。この現象（流れ）が観測されたのだと思います。この現象は3D では 3D でのシミュレーションの必要性があります。
＜電気泳動( v el )とな のシミュレーションの必要性があります。
みだしに DNA が引か ③ DNA がポアに停滞中のポア③ DNA がポアに停滞中のポアを見積もると、実験イルが一致します 5b 。

球対称となれば、点アの例ですが、この。実験結果は、V gate

V gate < V 溶 液

離れる(緑)

次の 6 秒間はポアか向に流れる異方性の流(緑 v osm )の和をの見積もり。

3

細胞の動力学への取り組み

3.1 実験：鶏 embryo の羽毛形成
次は鶏の embryo と発生の実験です。鶏は孵化するまえに皮膚に羽毛のパターンが形成します。先ずは背中
midline と呼ばれる背骨のラインに羽毛が等間隔に並び、それから左右対称に細密の六角形パターンができます[13]。つまり最初から生える場所と順序が決められています。現在、この場所のメモリを erase した人工皮膚がつくれます [15]。具体的には embryo の皮膚の二つの組織である上皮と間葉を、手で分けて、更に1細胞にばらしてから、それぞれの細胞を再集合させて、最後に上皮と間葉をくっつけます。midline のような genetic な成長の場
所と順番の記憶がない、『完全に一様な人工皮膚』ができます。そこで、まずは一次元の人工皮膚をつくり、羽毛の成長を観測しました。そして特徴的な成長場所が一次元の境界近くであることを発見しました [16]。その結果を COMSOL でシミュレーションをしてみました。
3.2

考えています。実験(赤)性のある運動のシミュレーションをするために3D 環境
で COMSOL を走らせております。この執筆には間に合Vol.22, No.4 2017
いませんでしたが、今後をご期待ください。

27.94 s28.58 s28.80 s29.23 s

図7 COMSOL 予想と実験結果。ポアを大きく、DNA を長く
図７. COMSOL 予想と実験結果。ポアを大きく、たくさんポアに置くと、DNA の逆流が引き起こされる。
DNA を長くたくさんポアに置くと、DNA の逆流が実験結果はそれを再現した。
引き起こされる。実験結果はそれを再現した。

実験では DNA がポアに停滞することがあります、ポ実験では DNA がポアに停滞することがあります、ポアに詰まっているかもしれません。そんな状況で後かアに詰まっているかもしれません。そんな状況で後から DNAがポアに近づくと、すり抜けてポアに入る場合
ら DNA がポアに近づくと、すり抜けてポアに入る場合もあれば、先に停滞しているDNAに はじかれる よもあれば、先に停滞している DNAに“はじかれる”よ
うに観える場合もあります。数年前に
DNA
由来の電気
うに観える場合もあります。数年前に DNA由来の電気
9浸透流が実験的に観測されました [12]。ちなみにその論浸透流が実験的に観測されました。ちなみにその論
COMSOL による解析でした。そこで、私はによる解析でした。そこで、私は“はじく”
文も文もCOMSOLはじ
条件を、実験前に
COMSOL
で予測しました。ポア内にく 条件を、実験前に COMSOL で予測しました。ポア
DNA を仮想的に置いて、電気浸透流の流れを観て、DNA
内に DNA を仮想的に置いて、電気浸透流の流れを観が“はじかれる”ものでしょうか。
て、DNA が はじかれる ものでしょうか。

シミュレーション：Turing の反応拡散方程式による羽毛のパタンーン形成
BZ 反応や、生物の皮膚に繰り返しパターンの形成が起きます。これを Turing が、一様な系からパターンが形成するまでの時間発展として、拡散項を含む偏微分
方程式（反応拡散方程式）によりモデル化に成功しました [17]。このモデルを Gierer や Meinhard らが発生現象に特 化 さ せ て、 そ れ ぞ れ2変 数 に よ る Activator Inhibitor
model, Activator Substrate model などを提唱しました [18]。
前者は羽毛のパターン、後者は歯の並びのシミュレー
ションに成功しました。近年、変数である Activator やInhibitor に対応する分子の候補が見つかってきました。
シミュレーション上の拡散係数にマッチするほどの複数の細胞間を通り抜ける膜貫通の機構は不明な点が多く、興味深いところです [19]。
3.3

COMSOL では、ポアを大きく、DNA を“多く”ポアに進入させると、DNA 由来の電気浸透流で、他の DNA が計算工学
ポアに進入しない予想をしました（図７上）。そして実験で試しました。T4 という購入できる最長(56.4 μm)の DNA を、そしてポアの径を 200(??-5)

nm にしました。実
Vol.22, No.4 2017（22-5）

チュートリアル：Turing の反応拡散方程式のシミュレーション
それでは、Activator Inhibitor model と Activator Substrate
Vol.23, No.1 2018

COMSOL Multiphysics による計算科学工学 –生物系 (4)te model をシミュレーションしましょう。事

一様ではない安定解(パターン)となる係数のセCOMSOL Multiphysics
による計算科学工学 –生物系 (4)チュートリアルCOMSOL Multiphysics による計算科学工学 – 生物系（4）

二つのパラメータの拡散係数の比(例えばDA << –生物系 (4)COMSOL Multiphysics による計算科学工学
COMSOL Multiphysics による計算科学工学 –生物系 (4)
ate
model
をシミュレーションしましょう。事
それらの大きさに対応するドメインサイズを予
model をシミュレーションしましょう。事前に、一様でSubstrate model をシミュレーションしましょう。事がポアにて
一様ではない安定解(パターン)となる係数のセければいけません。数学的には、ここが一つのはない安定解（パターン）となる係数のセットや二つのbstratemodel
をシミュレーションしましょう。事前に、一様ではない安定解(パターン)となる係数のセ
観測されま
二つのパラメータの拡散係数の比(例えばDA <<DA<<DI）と、それらなります。定常解として、左辺を0 に、そしてパラメータの拡散係数の比（例えばに、一様ではない安定解(パターン)となる係数のセ
ットや二つのパラメータの拡散係数の比(例えばDAA <<
めて観測さ、それらの大きさに対応するドメインサイズを予を除く、A, I の増減に関わる係数のdet を計算
Dの大きさに対応するドメインサイズを予測しなければいII )と、それらの大きさに対応するドメインサイズを予トや二つのパラメータの拡散係数の比(例えばDA <<チュートリアルチュートリアル
ければいけません。数学的には、ここが一つのでのDNAの測しなければいけません。数学的には、ここが一つのけません。数学的には、ここが一つの学問となります。
くと、それぞれの大小関係がわかります。ここ)と、それらの大きさに対応するドメインサイズを予します。シミュレーションでは初期条件を図10上のよ
うにすると、kA=0（図10左カラム）、kA ≠0（図10右カラム）と綺麗な六角形の spot と縦ラインの stripe が現れます。シマウマですね。
COMSOL
Multiphysics
COMSOL
MultiphysicsによCOMSOL は、ナノポアで慣れ親しんだ希釈種輸送を
チュートリアルCOMSOL Mul
す。COMSOL
学問となります。定常解として、左辺を0に、そして定常解として、左辺を0に、そして拡散項を除く、A、I
選びました。8.2の Nernst-Planck equation から電気泳動なります。定常解として、左辺を0 に、そしてターンができればOK として、参照文献のパラメ
しなければいけません。数学的には、ここが一つのDNA の異方拡散項を除く、A, I の増減に関わる係数の det を計算る過程においてCOMSOLのセット連続して計算しました。計測時間は5 分程度で る過程において
COMSOLの経
セット連続して計算しました。計測時間は5 分程度で
14, 16の増減に関わる係数のdetを計算していくと、それぞれと対流の項を外して、 R のリアクション項に、上記の
を除く、
A,
I
の増減に関わる係数のdetを計算
そのまま使用します。
問となります。定常解として、左辺を0 に、そしてる過程において C
セット連続して計算しました。計測時間は 5 分程度でめに３D 環していくと、それぞれの大小関係がわかります。ここくパラメトリックスイー
した。初期条件はCOMSOL
既存のRandom関数を一様なくパラメトリックスイープをした。初期条件はCOMSOL
既存のRandom関数を一様なkA をパラメトリッの大小関係がわかります。ここではパターンができれば右辺の二項目以降をかきます。またくと、それぞれの大小関係がわかります。ここ
散項を除く、
A,
I の増減に関わる係数のdet を計算
ctivator
Inhibitor
model
くパラメトリック
した。初期条件は COMSOL 既存の Random 関数を一様なには間に合ではパターンができればOK として、
参照文献のパラメ
し、このパラメトリック
微小ノイズに、図１０の空間異方性は、上記のノイズ
OK として、参照文献のパラメータをそのまま使用しま
クスイープで3セット連続して計算しました。計測時間
し、このパラメトリックスイ
微小ノイズに、図１０の空間異方性は、上記のノイズ
14,
14, 1616ていくと、それぞれの大小関係がわかります。ここターンができればOK として、
参照文献のパラメ
ator
を A, Inhibitor
を I として、。
ータをそのまま使用します。
し、このパラメト
微小ノイズに、図１０の空間異方性は、上記のノイズ
す [18, 20]。14,は5分程度でした。初期条件はCOMSOL 既存のRandom
ラメータを関数として扱
exp(cx 2 )をかけて初期条件としました。生命のパラメータを関数として扱えばに にexp(cx 2 )をかけて初期条件としました。生命のパ①
Activator
model
はパターンができれば
OK 16
として、
参照文献のパラメ
そのまま使用します。 Inhibitor
ラメータを関数とに exp(- cx 2 )をかけて初期条件としました。生命のパ
14,
16 Inhibitor①Activator
Activator
Inhibitor
model
関数を一様な微小ノイズに、図10の空間異方性は、上を A,を I として、各パラメータを同時にターンは拡散反応方程式に準ずるようですが、秩序あ各パラメータを同時に、最ターンは拡散反応方程式に準ずるようですが、秩序あ
タをそのまま使用します。
各パラメータを
ターンは拡散反応方程式に準ずるようですが、秩序あActivator
Inhibitor
model
2
A、Inhibitor を I として、
 A をAmodel
Activator
記のノイズにexp（-cx2）をかけて初期条件としました。
2
きると、最近知りました
るパターンを形成する場合には、初期条件か、境界条
Activator
Inhibitor
きると、最近知りました。Cるパターンを形成する場合には、初期条件か、境界条
きると、最近知りるパターンを形成する場合には、初期条件か、境界条
vator AA  A2
D A を A,
A Inhibitor2 を I として、孵化するま生命のパターンは拡散反応方程式に準ずるようですAkAA A Iを
A
Activator を A,
Inhibitor
レンドリーです。
件（ドメインの形）に空間的な異方性が組み込まれてレンドリーです。
件（ドメインの形）に空間的な異方性が組み込まれて
1
レンドリーです。
2 I として、先ずは背中
 D A A   A A   A 件（ドメインの形）に空間的な異方性が組み込まれて2が、秩序あるパターンを形成する場合には、初期条件
いるように思えます。
いるように思えます。
t1 kA A I
いるように思えます。
等間隔に並
2
2
2
か、境界条件（ドメインの形）に空間的な異方性が組み
nAーンができ
A
2
D

I


A


I



2
A
n n

I
I  A I
I 2 A
I  2A

AD
2


nA
＋
A








D
I
A
I
A




A
A
A
込まれているように思えます。
2
II  A I  I
 D  A
n
t
が決められ

が決められ At
se した人工1
t k A AI

1  k A A2 I

A

図８ A と I

皮膚の二つ
I I 2 2 I 2 2 I
I DI  AI,I 
更に１細胞以外のパラメータは定数
I A  I I  I
外のパラメータは定数

t
させて、最
で、先に結果を報告します。

 D  I  A  I 

Aに結果を報告します。

あることををシャープにします。
図９は、８に拡散係数D A <<
D I の効
拡散係数 D A <<D I の効によりピークのエッジ
ュレーショ
一
様
な
微
小
ノ
イ
ズ
(small
として、
A より
broad
て、 A より
broad
なな
I I
ープにします。
図９は、
fluctuation)からのパターン

fluctuation)からのパターン
ピークが確認できます。こクが確認できます。こ微 小 ノ イ ズ (small形成です。式の k AA を 0(図９左カ ラ ム )か ら 少 し ず つ 大 き くtion)からのパターン

違いによりピークのエッジ
拡散方程式
によりピークのエッジ

Activator Substrate model

Activator
Substratemodel
図８ A ②と I ② Activator
Substrate
model

図８
IとI I
図8
e のような
図８に拡散係数 D AA << D II の効
図８AAと
Aと
拡散係数
D I の効 A より broad な I、『完全に一D A <<果として、外のパラメータは定数
I 以外のパラメータは定数て、A より broad
な I
は一次元の
のピークが確認できます。こした。そしに結果を報告します。
の違いによりピークのエッジ、先に結果を報告します。
クが確認できます。こ

②

t t

②

Activator Substrate model

A2

 A2 SS

2S A SA
2
 AA  A A 2DSAAS A 
A A
2
2A
 DAD A2 
A A tA A 2  2A 
A1AAkAAA
kA A
t t
1 1kA A
2

S2 SD
AI2S S  1  k A 2    S
S S
2 2 St A
DI 
 DI    S A
S S 
k A2 A 2 St t
1 1kA A
ActivatorとSubstrate をそれぞれ   



A
n
A A
2 2 DcA1
D D A  A
At A 

 

t t
c3 c
I
 DI 
I I
2
I t
D2I 
Ic4 nA
c4 n
 DI t tn x , y ,t 0.0n xn,yx,,t y,t0 00.005

次は発生の実験

最先端の研究です
次は発生の実験＋モデ
次は発生の実験＋モデルの
図１１. A と S
分野では、シミュ

図１１.
Activator
図 ９ 時 間 (上 下 )。
図１１.
A とA Sと S
Activatorと とA, S とします。このモ図11 A と S
左から k A = 0, 0.1,
Substrate をそれぞれ最先端の研究です。こち最先端の研究です。こちらのをしている発生の分野では、シミュレーシ
分野では、シミュレーション

プは少ないのです
シャープにします。
図９は、Substrate をそれぞれするとspot-stripe の相変化をしている発生の実験グープにします。
図９は、をしている発生の実験グルー
デルでは、 A と S が負0.4。
す。
式の k A を 0(図９左率は 8 割を越えてンの形成が様 な 微 小 ノ (図９、中、右カラム)が観ら
イ ズ (small 図 ９ 時 間 (上 下 )。
S とします。このモ
A, SA,とします。このモの相関となります、(図プは少ないのですが、CO
なか微
小ノず
イつズ大(small
プは少ないのですが、COMSO
ら
少
し
き
く
の Michon の提唱
パターンが左から k Aembryo= 0, 0.1,
れます。さて、実際の鶏のuctuation)からのパターン
デルでは、
S が負
デルでは、
A とA Sとが負
１１,１２)。ピークに率は8割を越えています
ation)からのパターン
率は8割を越えています。こspot-stripeの相変化
デ ル は 、 羽 毛含む偏微分成です。式の k Aでは、上記で説明したを 0(図９左 0.4。midline にの相関となります、の相関となります、(図(図動きもあり、隙間をみ
k
図9
時間
（上下）
。左から
=0、0.1、0.4。
A
図
９
時
間
(上
下
)。
の
Michonの提唱したこ
す。
式の
k
を
0(図９左
成功しまし
の
Michonの提唱したこのモ
hexagonalのパタ
A
中、右カラム)が観ら
ラ ム )か ら 少 し一列に羽毛が並んでから、左右対ず つ 大 き く 図 ９ 時 間 (上 下 )。
１１,１２)。ピークにつけては、新たにピー
１１,１２)。ピークに
左から
k
=
0,
0.1,
A
が発生現象
称的に生成します。シミュレーシ生える様子を、上
デは
ル、
は羽
、
羽が
毛完
が璧
完な
璧
ム
)か ら 少 し ず つ 大 き
くデ ル毛。さて、実際の鶏のembryoると spot-stripe の相変化
左からクが生成するために密0.4。k A = 0, 0.1,
動きもあり、隙間をみ
Activator
動きもあり、隙間をみョンでは初期条件を図１０上のよ2種類の細胞を介
hexagonalのパターンとspot-stripeの相変化
hexagonalのパターンとして上記で説明した
midline に 0.4。
図９、中、右カラム)が観ら
度が A-I より一様です。
odel などをうにすると、 k AA = 0(図１０左カラ
つけては、新たにピー
つけては、新たにピーをそれぞれ成長生える様子を、上皮、間、中、右カラム)が観ら生える様子を、上皮、間葉の羽毛が並んでから、左右対ます。さて、実際の鶏のembryo
後者は歯のム)、 k AA ≠ 0(図１０右カラム)とドメインサイズが時間
クが生成するために密BMP7 と BMP2 を候クが生成するために密2 種類の細胞を介して、。さて、実際の鶏のembryospot
2 種類の細胞を介して、
A とA
近年、変数
生成します。シミュレーシ
は、上記で説明した
midlineに と縦ライン綺麗な六角形の的に大きくなる系(成立てて式を立てま度がA-Iより一様です。
度がA-I
より一様です。
分子の候補の stripeが現れます。
シマウマで
をれ
そぞ
れれ
ぞ成
れ長
成に
長関
にわ
関る
わ
上記で説明した
midlineにをそ列に羽毛が並んでから、左右対は初期条件を図１０上のよ長)ではピークの数が
図１２. 時間(上下)。左
れまでの概念を覆ドメインサイズが時間上の拡散係
ドメインサイズが時間
すね。
から
k
=
0,
0.1,
0.4。
BMP7
と
BMP2
を候補分子
的に生成します。シミュレーシ
羽毛が並んでから、左右対
A
BMP7
と
BMP2
を候補分子と見密度を一定に保つようると、 k A = 0(図１０左カラ舞いをシミュレー
抜ける膜貫
図１０. 初条大COMSOL は、ナノポアで慣れ親し
的
にき
大く
きな
くる
な系
るに増えるので、口の大
系 (成
的期に
(成
17 17
ンでは初期条件を図１０上のよ立てて式を立てました
生成します。シミュレーシ
立てて式を立てました。こ≠
0(図１０右カラム)と15
15
件を
midlineに。
験結果で確かめま
図12 時間（上下）。左から kA=0、0.1、0.4。
ろです 。
んだ希釈種輸送を選びました。８．
長
)
で
は
ピ
ー
ク
の
数
が
図１２.
時間(上下)。左
長
)
で
は
ピ
ー
ク
の
数
が
図１２.
時間(上下)。左きさと歯の数の増加などに、このモデルが用いられま
にすると、
k
=
0(図１０左カラ
k
=
0
と
0.4。
A
A
れまでの概念を覆す
BMP
は初期条件を図１０上のよAれまでの概念を覆す
BMPの振六角形の spot２のと縦ライン
成長を、細胞の密Nernst-Planck equation からから= 0.1,
0, 0.1,0.4。
16から
k A =k A0,
0.4。
A、S とします。この
密度を一定に保つよう
ActivatorとSubstrateをそれぞれ
密度を一定に保つよう
す
。
)、
k
≠
0(図１０右カラム)とA
舞いをシミュレーション
ると、
k A = 0(図１０左カラ散方程式の電シマウマで
気泳動と対流の項を外し
舞いをシミュレーションと実ます(論文では二peが現れます。
モデルでは、Aと S が負の相関となります、（図11、12）。
に増えるので、口の大これらのモデルは
2 変数によるシミュレーションで、に増えるので、口の大
spot
と縦ライン
て、
”R”のリアクション項に、上記の右辺の二項目以
ミュレーションす験結果で確かめました。
k麗な六角形の験結果で確かめました。このA ≠ 0(図１０右カラム)とピークに動きもあり、隙間をみつけては、新たにピーク
手計算にて一様ではない安定解(パターン)となる各定
Activator
降をかきます。また
k AA をパラメトリックスイープで３きさと歯の数の増加などに、このモデルが用いられま
stripeが現れます。
シマウマできさと歯の数の増加などに、このモデルが用いられましましょう。図１
図10
初期条件を
midlineに。k=0と0.4。
成長を、細胞の密度としspot と縦ライン
図１０. 初期条 A
成長を、細胞の密度として3
L六角形のは、ナノポアで慣れ親し16が生成するために密度がA-I より一様です。ドメインサ
数のパラメータの大小比を予測できますが、変数が
3
16
ね。
す
。
す
。
midline から順番
件をmidlineに。
ます(論文では二種類の細ipeが現れます。
シマウマでます(論文では二種類の細胞
釈種輸送を選びました。８．
つ以上になると、本が一冊書けます。余談ですが、私
イズが時間的に大きくなる系（成長）ではピークの数が 10b 。 A の 勾 配 に. 初Vol.22,期 条 No.4 2017
COMSOL は、ナノポアで慣れ親し(??-6) A、I 以外のパラメータは定数で、先に結果を報告しこれらのモデルは
2 変数によるシミュレーションで、これらのモデルは
2 変数によるシミュレーションで、
k A =図0１と０0.4。
ミュレーションする場合
ミュレーションする場合は原nst-Planck equationから件を
midlineに。
は COMSOL の力技で 3 変数のモデルつくった経験があり密度を一定に保つように増えるので、口の大きさと歯のま す。 図8に 拡
散係数
DA<<D
果とし
て、A
よ
り
Iの効
だ希釈種輸送を選びました。
８．
chemotaxis)と、そ手計算にて一様ではない安定解(パターン)となる各定
手計算にて一様ではない安定解(パターン)となる各定
図
１
０
.
初
期
条
[20]
しし
まょ
しう
ょ。
う図
。１
図３
１の
３結
の果
結
OL
は、ナノポアで慣れ親し
k
=
0
と
0.4。
し
ま
A
動Nernst-Planckと 対 流 のbroad
項equationをしからます。 A と I と新しく M を加えました。実験的に、三数の増加などに、このモデルが用いられます。
I のピークが確認できます。この違いによりな外の献が、結果として数のパラメータの大小比を予測できますが、変数が3数のパラメータの大小比を予測できますが、変数が3
件をmidlineに。
midlineから順番に出来釈種輸送を選びました。
８．
midlineから順番に出来るこつ分子によるフィードバックを予想して、後はパラメ”のリアクション項に、これらのモデルは2変数によるシミュレーションで、気 泳 動 と 対 ピークのエッジをシャープにします。図9は、一様な微流 の 項 を上記の右辺の二項目以外し秩序化を促進しまk A = 0 と 0.4。
つ以上になると、本が一冊書けます。余談ですが、私
つ以上になると、本が一冊書けます。余談ですが、私
10b 10b 。
トリックスイープ
allcombinationsの力任せです。
そ
A
の
勾と
配
に
よ動
りこ
動
ernst-Planck
equation
から
。
A
の
勾
に
く
手計算にて一様ではない安定解（パターン）となる各定 配
きます。またk小ノイズ
A をパラメトリックスイープで３（small ﬂuctuation）からのパターン形成です。
、”R”のリアクション項に、上記の右辺の二項目以
こよ
ろり
で
、
の
は COMSOLの力技で3 変数のモデルつくった経験がありはCOMSOLの力技で3変数のモデルつくった経験がありして実験結果と共に投稿しましたが、見事に力不足とchemotaxis)と、その自己泳
動 と 対 流 式のの項
外
し数のパラメータの大小比を予測できますが、変数が3つ
chemotaxis)と、その自己増殖を0（図9左カラム）から少しずつ大きくするとCOMSOL を使ってAををかきます。また
k Akをパラメトリックスイープで３ます。
I と新しくM を加えました。実験的に、三ます。
A とA Iとと新しく
M を加えました。実験的に、三判断が査読にて下りました。
「三つ目の M の式は A献が、結果としてピークとI以上になると、本が一冊書けます。余談ですが、私はR”のリアクション項に、献が、結果としてピークの動
spot-stripe の 上記の右辺の二項目以
相 変 化（図9、 中、 右 カ ラ ム）が 観 ら れ ま
COMSOL を用いておVol.22, No.4 2017
つ分子によるフィードバックを予想して、後はパラメへ影響を与えない。それでいて
A と I から M 用のパタ
つ分子によるフィードバックを予想して、後はパラメ
COMSOLの力技で3変数のモデルつくった経験がありりませんが、パラ
embryo では、上記で説明した
秩序化を促進しました。
きます。また す。さて、実際の鶏のk A をパラメトリックスイープで３秩序化を促進しました。
Vol.22, No.4 2017
ーンを。それは意味がない」と。さすが専門家です。
トリックスイープ
allcombinations
トリックスイープ
allcombinationsの力任せです。
そ そ
Mの力任せです。
ます。Aと I と新しくを加えました。実験的に、三つ

の
n
性の
midline に一列に羽毛が並んでから、左右対称的に生成
とろ
こで
ろ、
で
、の
こ図
図
１A
３
とこ
こ
１
３
は
まさにその通りでした。しかし、私は、モデルを立てして実験結果と共に投稿しましたが、見事に力不足として実験結果と共に投稿しましたが、見事に力不足と期条件は関数で定COMSOLを使って計算した
COMSOLを使って計算したもVol.22, No.4 2017
計算工学
（22-6）
Vol.23,
2018
判断が査読にて下りました。
「三つ目の
M の式は
I No.1
判断が査読にて下りました。
「三つ目の
M の式は
A とA IとCOMSOLを用いており、論
COMSOLを用いており、論文にへ影響を与えない。それでいて
I から
M 用のパタ
計算工学(??-7)
へ影響を与えない。それでいて
A とA Iとから
M 用のパタりませんが、
パラメータはりませんが、
パラメータはすーンを。それは意味がない」と。さすが専門家です。
ーンを。それは意味がない」と。さすが専門家です。
n
A は
性の
 n
A は性の
PDEPDEを使まさにその通りでした。しかし、私は、モデルを立てまさにその通りでした。しかし、私は、モデルを立て期条件は関数で定義しま
期条件は関数で定義しまして   

チュートリアルCOMSOL Multiphysics による計算科学工学 – 生物系（4）分子によるフィードバックを予想して、後はパラメト
リックスイープ all combinations の力任せです。そして実験結果と共に投稿しましたが、見事に力不足と判断
COMSOL Multiphysics
による計算科学工学–生物系 (4)
M の式は A と I へ影が査読にて下りました。
「三つ目の響を与えない。それでいて A と I から M 用のパターン
る過程において COMSOL の経験値を上げました。効率よを。それは意味がない」と。さすが専門家です。まさにくパラメトリックスイープをする方法などです。ただCOMSOL Multiphysics による計算科学工学 –生物系 (4)その通りでした。しかし、私は、モデルを立てる過程
し、このパラメトリックスイープにおいて、複数のパにおいて COMSOL の経験値を上げました。効率よくパる過程において COMSOL の経験値を上げました。効率よは 5 分程度でラメータを関数として扱えば、つまり、
ラメトリックスイープをする方法などです。ただし、
くパラメトリックスイープをする方法などです。ただ関数を一様な
各パラメータを同時に、最適に変化させることもで
このパラメトリックスイープにおいて、複数のパラ
し、このパラメトリックスイープにおいて、複数のパ
上記のノイズ
きると、最近知りました。COMSOL
は、とても USER フ
メータを関数として扱えば、つまり、ラメータを関数として扱えば、つまり、た。生命のパレンドリーです。
各パラメータを同時に、最適に変化させることもで各パラメータを同時に、最適に変化させることもですが、秩序あきると、最近知りました。COMSOL
は、とてもは、とてもUSER フ件か、境界条きると、最近知りました。COMSOL
USER

k
n
N

n
t

t

レンドリーです。
組み込まれてn
フレンドリーです。
n







   nA＋
- k n nt kt nN  n t  t tn   nA2＋ n- k n nt  t t c n 1  c 2 A A D A 2 A  1kAAA c3  A22 1 c1In1  c 2 A 2 t D  A
k A
t

A







A



c  A 1  I 
2



拍動を単

機能は、

個々の拍

A

観測でき

体は、部
図14

最先端の研究です。こちらの最先端の研究です。こちらの１. A と S

分野では、シミュレーション
分野では、シミュレーションをしている発生の実験グルーをしている発生の実験グルー

プは少ないのですが、COMSOL

プは少ないのですが、COMSOL
率は 8 割を越えています。こ率は 8 割を越えています。この Michon の提唱したこのモ

デ ル は 、 羽 毛 が 完 璧 なの Michon の提唱したこのモ

hexagonal
デ ル は 、 羽 毛が 完 璧のパターンとしてな生える様子を、上皮、間葉のhexagonal のパターンとして2 種類の細胞を介して、A と I生える様子を、上皮、間葉のをそれぞれ成長に関わる

2 種類の細胞を介して、A とを候補分子と見I
BMP7 と BMP2
図 １ ３ . A(右 )と ｎ

れまでの概念を覆す BMP の振BMP7 と BMP2 を候補分子と見(左)0.1, 0.4。
舞いをシミュレーションと実17立てて式を立てました。こ
図13 A（右）と n（左）図 １ ３ . A(右 )と ｎ
験結果で確かめました。このモデルは最終的な羽毛のれまでの概念を覆す
BMP の振が用いられま(左) 3 つ目の方程式も用いてい
成長を、細胞の密度として次は発生の実験＋モデルの最先端の研究です。こち
舞いをシミュレーションと実ます(論文では二種類の細胞として4 つの方程式)。シらの分野では、シミュレーションをしている発生の実
レーションで、験結果で確かめました。このモデルは最終的な羽毛のミュレーションする場合は原著のパラメータを参考に
験グループは少ないのですが、COMSOL 率は8割を越ン)となる各定
し ま し ょ う 。3図つ目の方程式も用いてい
１３の結果は、六角形のパターンが成長を、細胞の密度としてえています。この Michon の提唱したこのモデルは、羽が、変数が 3
midlineから順番に出来ることを完璧に再現しましたます(論文では二種類の細胞として 4 つの方程式)。シ
談ですが、私10b 毛が完璧な hexagonal のパターンとして生える様子を、。 A の 勾 配 に よ り 動 く 細 胞 n( 細 胞 の 走 化 性 ,ミュレーションする場合は原著のパラメータを参考にた経験があり上皮、間葉の2種類の細胞を介して、Aと I をそれぞれchemotaxis)と、その自己増殖、そしてA の増加への貢しましょう。図１３の結果は、六角形のパターンが実験的に、三成長に関わる BMP7と BMP2を候補分子と見立てて式を献が、結果としてピークの動きを抑制して、空間的な、後はパラメmidlineから順番に出来ることを完璧に再現しました
秩序化を促進しました。
立てました [21]。これまでの概念を覆す BMP の振舞いを0b任せです。
。 A のそ勾 配 にと
よこりろ動
くこ
細の胞
(細
胞私のが走
性式,を 参 考 に
で、
図 １n３
は、
論化
文の
シミュレーションと実験結果で確かめました。このモ事に力不足とCOMSOLを使って計算したものです。この論文自体もA の増加への貢
chemotaxis)と、その自己増殖、そしてデルは最終的な羽毛の成長を、細胞の密度として3つ目の式は A と I

COMSOL を用いており、論文にフィジクス等の説明はあ献が、結果としてピークの動きを抑制して、空間的なの方程式も用いています（論文では二種類の細胞として
ら M 用のパタ
りませんが、パラメータはすべて SI にあります。走化
秩序化を促進しました。
4つの方程式）。シミュレーションする場合は原著のパ専門家です。
性の   nA は PDE を使い簡単にかけました。初モデルを立てと こ ろ で 、 こ期条件は関数で定義しまして、時間による関数の変化のラメータを参考にしましょう。図13の結果は、六角形
図１３は、私が論文の式を参考に



COMSOL を使って計算したものです。この論文自体ものパターンが midline から順番に出来ることを完璧に再[14]。A の勾配により動く細胞（細胞の走化COMSOL を用いており、論文にフィジクス等の説明はあ(??-7) 現しました
Vol.22, No.4 n
2017

りませんが、パラメータはすべて SI にあります。走化

心筋細

ーレ上再

次は発生の実験＋モデルの次は発生の実験＋モデルの



て物理的

ます。こ



17
を そ れ ぞ れ 成 立てて式を立てました
長に関わる
。こ

４．心筋

４．1 実

細胞の電

3
I
 D I  2 I Ic 4 nA 2 2k I I
 D I  I  c 4 nA 2  k I Ittn x , y ,t  0 n0x.,005
 200
x 2  200 x 2 
y ,t 02 exp
 0.005
 2exp間(上下)。左性、chemotaxis）と、その自己増殖、そして A の増加への貢献が、結果としてピークの動きを抑制して、空間的な秩序化を促進しました。
と こ ろ で、 こ の 図13は、 私 が 論 文 の 式 を 参 考 にCOMSOL を使って計算したものです。この論文自体もCOMSOL を用いており、論文にフィジクス等の説明はチュートリアルありませんが、パラメータはすべて SI にあります。走化性の∇・(χn∇A) は PDE を使い簡単にかけました。初は step( t )を用いました。
期条件は関数で定義しまして、時間による関数の変化
最後に私たちの実験です。モデルの詳細は投稿中のは step(t) を用いました。
ため省略しますが、実験結果はシンプルで、どちらか
最後に私たちの実験です。モデルの詳細は投稿中のというとシミュレーションの結果を、実験が忠実に再ため省略しますが、実験結果はシンプルで、どちらか現しました。
というとシミュレーションの結果を、実験が忠実に再現しました。

羽毛の発生の断面図とその模式図（左図）。Turing モデ

図１４.
羽毛の発生の断面図とその模式図(左図)。
ルを立てて幅を変え。二つの変数にて、境界条件を、
Turing 上：二つとも流束0、下：ノイマン＋ディリクレのモデルを立てて幅を変え。二つの変数にて、mix（下）。
境界条件を、上：二つとも流束
0、下：ノイマン＋
ディリクレの mix (下)。
一般的に反応拡散モデルの境界条件は粒子保存の流束0（ノイマン境界条件）です。しかし、1D 的に細長い
2D一ドメインでシュミレーションすると、当然ながら般的に反応拡散モデルの境界条件は粒子保存のエッジにピークができます。実験では真ん中一列にし流束
0(ノイマン境界条件)です。しかし、1D 的に細長か羽毛が生えず、図14の断面図では、皮膚のエッジ（左い 2D ドメインでシュミレーションすると、当然ながら右の端）から一定の間隔に羽毛が生えます。羽毛の間隔
エッジにピークができます。実験では真ん中一列にしは A-I の場合、I の拡散係数で決まります。そこで、I をか羽毛が生えず、図１４の断面図では、皮膚のエッジ
ノイマン境界条件で、A はエッジにピークが出来ない(左右の端)から一定の間隔に羽毛が生えます。羽毛のように流束≠0、あるいはディリクレ境界条件の A=0
間隔は A-I の場合、I の拡散係数で決まります。そこと、境界条件を mix にするとシミュレーションは実験で、I をノイマン境界条件で、A はエッジにピークが出結果を再現しました。図14のような、実験家の案らし来ないように流束≠0、あるいはディリクレ境界条件のいドメインの形状や、mix の境界条件などは COMSOL
Aを利用すれば、直感的にシミュレーションを構築でき= 0 と、境界条件を mix にするとシミュレーション
は実験結果を再現しました。図１４のような、実験家
走らせることができます。
の案らしいドメインの形状や、mix の境界条件などは
COMSOLを利用すれば、直感的にシミュレーションを構
3.4 実験結果との考察：人工皮膚上の羽毛と反応拡散のシミュレーションから予測できること
築でき走らせることができます。
二つの A-I と A-S モデルは、生物の体表に現れる模様
や毛の配置のパターン形成をシミュレーションとして３．
４ 実験結果との考察：人工皮膚上の羽毛と反応拡再現しました。しかし、A-I
では二つのピークが合体し散のシミュレーションから予測できることたり、A-S
では、ピークが動き回ったりと、観測されな二つの A-I と A-S モデルは、生物の体表に現れる模い現象も見受けられます。しかし、三つ目の式（細胞）様や毛の配置のパターン形成をシミュレーションとしにより、実際に観測される羽毛の生成過程の忠実な再て再現しました。しかし、A-I では二つのピークが合現に成功しました。これによりシミュレーションの各体したり、A-S では、ピークが動き回ったりと、観測パラメータから、各分子の細胞間の移動について新たされない現象も見受けられます。しかし、三つ目の式な見識が加わりました。我々の実験では、発生の初期(細胞)により、実際に観測される羽毛の生成過程の忠条件を一様にすることに成功し、反応拡散に忠実に羽実な再現に成功しました。これによりシミュレーショ
毛が発生することもわかりました。生命科学の飛躍に
ンの各パラメータから、各分子の細胞間の移動につい
大きくシミュレーション（COMSOL）が貢献しています。
て新たな見識が加わりました。我々の実験では、発生

異方性が

の伝搬は

みます。
す。

私は、

のシャー

もつ系へ

すすめて

ら与えた
４．２

シミュ

を模した

タートす

C
I
I
I

デル化し

チャネル

明します

す。細胞

それぞれ

に変化す

h は時間

で開閉す

このよう

計算工学
性の   nA は
PDE を使い簡単にかけました。初（22-7）
Vol.23, No.1 2018の初期条件を一様にすることに成功し、反応拡散に忠は、個々

期条件は関数で定義しまして、時間による関数の変化

実に羽毛が発生することもわかりました。生命科学のるため、飛躍に大きくシミュレーション(COMSOL)が貢献してい

ャネルを

ます。

純な実験

Vol.22, No.4 2017

チュートリアルCOMSOL Multiphysics によ

要があります。 V を V ( t , x ,COMSOL Multiphysics による計算科学工学 –生物系 (4)

投稿中の

どちらか

忠実に再

V
を第一式の右辺に加えます。
 V VCOMSOL  1Multiphysics
 V   w による計算科学工学
 Iap
チュートリアル– 生物系（4）４．心筋細胞への数値的アプローチ
でかけません。そしてドメイt
COMSOL Mu
４．1 実験: 心筋細胞の拍動とその空間的伝搬、そしません。
wチュートリアル  V  w 
4 心筋細胞への数値的アプローチて物理的刺激の励起する拍動について
t4.1
実験：心筋細胞の拍動とその空間的伝搬、そして要があります。 V
心筋細胞はexcitable な細胞として周期的な自律
４．３ チュートリアル：Fit
V
物理的刺激の励起する拍動について
拍動を単一細胞にばらしても行います。心臓のポンプを第一式の右辺に V V   1  V   w  Iap空間によるシミュレーション
一 方 、 FitzHugh-Nagumo
機能は、この心筋細胞と、自律拍動をしないが、心筋心筋細胞は excitable な細胞として周期的な自律拍動でかけません。そtmodel は 、 パ ル ス 状 に細胞の電気的な結合を補助する繊維芽細胞が担っていを単一細胞にばらしても行います。心臓のポンプ機能ません。
私は COMSOL
希釈種輸送をwます。これらの細胞を、単一細胞にバラして２D
シャexcite する細胞の電位と、  V   w 
は、この心筋細胞と、自律拍動をしないが、心筋細胞使います。先ずは、定義でドt
ーレ上再結合させると、一様な系として細胞集団の その電位と相反する、たっの電気的な結合を補助する繊維芽細胞が担っていま

４．３ チュート
メインポイントプローブを置
空間によるシミュいておきます。リアルタイム

個々の拍動と、時間経過と共に起こる拍動の同期化をた二つの微分方程式から成す。これらの細胞を、単一細胞にバラして2D シャーレ

21 FitzHugh-Nagumo
一方、
観測できます。まさに蔵本モデルです。一方、心臓本り、とても単純です。私

上再結合させると、一様な系として細胞集団の個々ので拍動の周期的振動を確認で拍動と、時間経過と共に起こる拍動の同期化を観測で
私は COMSOL 希
きます(図１７)。ジオメトリ
異方性が組み込まれています。例えば、拍動シグナルでコード化するほうが早い
excite する細胞の電位と、きます。まさに蔵本モデルです。一方、心臓本体は、使います。先ずはでは、円のドメインを 4 つにの伝搬は、ペースメーカー細胞から、特定の方向に進
その電位と相反する、たっくらいです。位相空間の振部位ごとに物理特性が異なり、細胞集団による異方性
メインポイントプ
わけます (図１６右下)。
初期
みます。もし、その方向性が崩れると不整脈になりまた二つの微分方程式から成動による軌跡は、それぞれが組み込まれています。例えば、拍動シグナルの伝搬
いておきます。リ
条件に空間異方性を与えるた
図15 V-w
す。
21位相と V、w の時間変化。
図。私
１ ５ . V-w 位 相に式の左辺をり、とても単純です
0 として、図は、ペースメーカー細胞から、特定の方向に進みまで拍動の周期的振めです。また小さい円はペー
私は、この心臓の系(心筋細胞＋繊維芽細胞)を 2Dと V, w の時間変化。
の場合、COMSOL より matlabする細胞の電位と、その電位と相反する、たった二つ
１５のように、Vの３次関す。もし、その方向性が崩れると不整脈になります。
きます(図１７)。
スメーカー用です。例えのシャーレ上で、一様(homogeneous)な系から異方性を数と線形のwでコード化するほうが早いとの間を振動の微分方程式から成り、とても単純です [25]。私の場 では、円のドメイ
私は、この心臓の系（心筋細胞＋繊維芽細胞）を2Dのもつ系への変化を実験的再現することを目的に研究をば I ap ＝I 0 δ(x, y)のようくらいです。位相空間の振する緑の曲線にて見積もることができます。速い時間
合、COMSOL より matlab でコード化するほうが早いく わけます (図１６
シャーレ上で、一様（homogeneous）な系から異方性をすすめています。例えばメカニカルな刺激を一方向かにとれば、拍動波面の湧動による軌跡は、それぞれ変化の V と遅い時間変化のw の大小により、右辺のらいです。位相空間の振動による軌跡は、それぞれにもつ系への変化を実験的再現することを目的に研究をら与えたら異方性は形成するのでしょうか。
条件に空間異方性
き出しを限定できます。
図 １ ５ . V-w
位相に式の左辺を0 として、図
正負が周期的に変わるからです。これらの式の各パラ式の左辺を0として、図15のように、Vの3次関数と線
４．２
シミュレーション：FitzHugh-Nagumo model
すすめています。例えばメカニカルな刺激を一方向かめです。また小さ
初期条件に位相差とし
V, w の時間変化。
１５のように、V の３次関 ε は と
ε<< 1形の w との間を振動する緑の曲線にて見積もることがシミュレーションでは、拍動シグナルの電位パルスメータの単位は“適当”ですが、例えばら与えたら異方性は形成するのでしょうか。
スメーカー用ですて図１６右下のようにす
数と線形のw との間を振動 V と遅い時間変化のw の大 ば I ＝I δ(x, yIapできます。速い時間変化のを模した微分方程式の FitzHugh-Nagumo model からスとして遅い時間変化(弱いカップリング)を与え、ap0
ると、半径の長さの波面する緑の曲線にて見積もることができます。速い時間
は拍動を励起する流入電流として、それぞれ生理学的小により、右辺の正負が周期的に変わるからです。こ にとれば、拍動波
タートするのですが、その前に細胞膜のチャネルをモ
4.2 シミュレーション：FitzHugh-Nagumo model
が時計回りに回転し、速変化のV と遅い時間変化のw の大小により、右辺のIap はシミュレーに意味があると考えられます。このれらの式の各パラメータの単位は適当 ですが、例え き出しを限定でき
さが r に依存しないこと
dV
正負が周期的に変わるからです。これらの式の各パラション時に、振動の条件を引き起こすパラメータとしばεはε <<1として遅い時間変化（弱いカップリング） 初期条件に位相C
I I I
I
から、螺旋の波面伝搬がK
Na
Cl
メータの単位は“適当”ですが、例えばε は ε<< 1dtてとても有益です。
を与え、Iap は拍動を励起する流入電流として、それぞて図１６右下のよ起こります(図１８右上)。
として遅い時間変化(弱いカップリング)を与え、Iap
さて、私の実験では、心筋細胞の初代培養を
2 次元
Iap
れ生理学的に意味があると考えられます。このはシ ると、半径の長さ
I  g n4 V  V
先ほどの I ap ＝I
0 δ(x, y)KKKは拍動を励起する流入電流として、それぞれ生理学的のシャーレ上で行い、細胞集合体の拍動と、その空間
ミュレーション時に、振動の条件を引き起こすパラ が時計回りに回転
は波面の湧き出しを中心に意味があると考えられます。このIap はシミュレー
的伝搬を観測しています。ゴールはその伝搬の機構とメータとしてとても有益です。
I g m 3h V  V
さが r に依存しなとして円状に波面が広がNa
NaNa
ション時に、振動の条件を引き起こすパラメータとし同期化、そして細胞のメカニカルな動きの拍動と、拍
さて、私の実験では、心筋細胞の初代培養を2次元のから、螺旋の波面ります。ペースメーカーI  g V  Vてとても有益です。
動伝搬に与える影響の理解です。つまり私の実験をシシャーレ上で行い、細胞集合体の拍動と、その空間的Cl
LCl起こります(図１ですね。COMSOL
ではδ関
さて、私の実験では、心筋細胞の初代培養を 2 次元
ミュレーションするには拍動を空間的に伝搬させる必
伝搬を観測しています。ゴールはその伝搬の機構と同先ほどのI ap ＝I 0数などの様々な関数に加シミュレーションでは、拍動シグナルの電位パルス
デル化した Hodgkin-Huxley model 18 を元に各種イオンのシャーレ上で行い、細胞集合体の拍動と、その空間
期化、そして細胞のメカニカルな動きの拍動と、拍動
は波面の湧き出し
えて、論理式も multiple でを模した微分方程式の FitzHugh-Nagumo model からス
チャネルの特性を含む Luo Rudy model 19 について説
的伝搬を観測しています。ゴールはその伝搬の機構と伝搬に与える影響の理解です。つまり私の実験をシとして円状に波面件を試せます。
Iap ≠0 とする
タートするのですが、その前に細胞膜のチャネルをモ明します。生物学的には、より理想的なモデルなので同期化、そして細胞のメカニカルな動きの拍動と、拍ミュレーションするには拍動を空間的に伝搬させる必[22]
ります。ペースメて、式内のパラメータを変えデル化した
Hodgkin-Huxley
modelを元に各種イオンす。細胞外の電位を 0 V とした際に、細胞内の電位が動伝搬に与える影響の理解です。つまり私の実験をシy) を第一 ですね。COMSOL で要があります。V を V(t, x, y) として D ∇2V(t, x,
うな、拍動伝搬の速さが異なチャネルの特性を含む Luo Rudy model[23] について説明それぞれのイオン種別チャネルの開閉により、時間的
ミュレーションするには拍動を空間的に伝搬させる必数などの様々な関につなげられて実験条件にあします。生物学的には、より理想的なモデルなのでに変化する振る舞いを予測します。
このモデルの n , m ,できます。 えて、論理式もとした際に、細胞内の電位がそh す。細胞外の電位を0Vは時間の関数であり、更に、細胞の内外での電位差
件を試せます。I
近年では、心臓がメカニカで開閉するチャネルもあり、V の関数でもあります。
れぞれのイオン種別チャネルの開閉により、時間的に体は、部位ごとに物理特性が異なり、細胞集団によるの場合、
model
は、パルス状にCOMSOL
より matlab

左図)。
数にて、
イマン＋

子保存の

的に細長

然ながら

一列にし

のエッジ

。羽毛の

す。そこ

ークが出

界条件の

ーション

、実験家

件などは

ョンを構

と反応拡

現れる模

ョンとし

ークが合

と、観測

つ目の式

過程の忠

レーショ

動につい

















は、発生

このように
Hodgkin-Huxley model や Luo Rudy mode 20 n、m、変化する振る舞いを予測します。このモデルの拡散に忠は、個々のチャネルの物理的振る舞いを正しく記述すh は時間の関数であり、更に、細胞の内外での電位差で命科学のるため、結果としてパラメータが多くなり、個々のチ
開閉するチャネルもあり、V の関数でもあります。こ

献していャネルをブロックするような試薬を用いない限り、単のように Hodgkin-Huxley model や Luo Rudy mode[24] は、
純な実験観測の結果とシミュレーションを比べること(V, w):初期条件
(1,0
(0,0

て、式内のパラメ
クとして、拍動を安定化する
うな、拍動伝搬のています。つまり心臓を弾性につなげられて実も心臓内の電流（電位変化）できます。
す。そこで、COMSOLの versa

近年では、心臓の力を借りて挑戦しました。
個々のチャネルの物理的振る舞いを正しく記述するた(0,0.1)ができません。
クとして、拍動をめ、結果としてパラメータが多くなり、個々のチャネ(V, w):初期条件ています。つまり
ルをブロックするような試薬を用いない限り、単純なチュートリアルCOMSOL
Multiphysics
による計算科学工学–生物系
チュートリアルCOMSOL
Multiphysics
による計算科学工学–生物系(4) (4)も心臓内の電流（図１６. 反応速度に FHN 式の左辺を。w の方は拡散実験観測の結果とシミュレーションを比べることがで(1,0(??-8)
Vol.22, No.4 2017 係数を 0 に。螺旋を描く波形の伝搬をシミュレーシす。そこで、COM
きません。
2
計算工学
要があります。
Vをy )として要があります。
VをV右下図のようにドメインを分けて、( tV,( tx,, xy,)として
D∇D∇
V ( t2 V,( tx, yx), y )(??-9)(0,0ョンする場合は、

V
V VV  11VVwwIap
Iap
V
t
w
 VV  ww
t(0,0.1)を第一式の右辺に加えます。この式は(私は)、Matlab
図のような初期条件を。
を第一式の右辺に加えます。この式は(私は)、Matlabの力を借りて挑戦

でかけません。そしてドメインも動かさなければいけでかけません。そしてドメインも動かさなければいけ
ません。
ません。

図１６.
反応速度に
式の左辺を。w
の方は拡散
図16
反応速度に
FHNFHN
式の左辺を。wの方は拡散係数を0係数をに。螺旋を描く波形の伝搬をシミュレーションする0 に。螺旋を描く波形の伝搬をシミュレーシ
４．３
model
の時の時
４．３チュートリアル：FitzHugh-Nagumo
チュートリアル：FitzHugh-Nagumo
model
計算工学
(??-9)
場合は、右下図のようにドメインを分けて、図のよョンする場合は、
右下図のようにドメインを分けて、一 方、FitzHugh-Nagumo model は、 パ ル ス 状 に
excite空間によるシミュレーション
うな初期条件を。
空間によるシミュレーション
図のような初期条件を。
一方方、、 FitzHugh-Nagumo
一
FitzHugh-Nagumo
model
は
、
model は
、パ
パル
ルス
ス状
状に
に
私は
COMSOL
計算工学（22-8）
Vol.23, No.1 2018
私は
COMSOL希釈種輸送を
希釈種輸送を
exciteする細胞の電位と、
excite する細胞の電位と、使います。先ずは、定義でド使います。先ずは、定義でドその電位と相反する、たっその電位と相反する、たっ
メインポイントプローブを置
メインポイントプローブを置た二つの微分方程式から成た二つの微分方程式から成いておきます。リアルタイム
いておきます。リアルタイム
り、とても単純です 21
り、とても単純です 21。私。私
で拍動の周期的振動を確認でで拍動の周期的振動を確認でがあります。
V を V ( t , xによる計算科学工学, y )として D∇ 2 V (–生物系
t , x , y(4))
COMSOL Multiphysics

チュートリアル

第一式の右辺に加えます。この式は(私は)、Matlab
要があります。
V を V ( t , x , y )として D∇ 2 V ( t , x , y )かけません。そしてドメインも動かさなければいけ

チュートリアルを第一式の右辺に加えます。この式は(私は)、Matlab
せん。
でかけません。そしてドメインも動かさなければいけ
式の右辺に加えます。この式は（私は）、Matlab でかけません。
．３ チュートリアル：FitzHugh-Nagumo model の時ません。そしてドメインも動かさなければいけません。
間によるシミュレーション
４．３ チュートリアル：FitzHugh-Nagumo model の時
4.3 チュートリアル：FitzHugh-Nagumo model の時
空間によるシミュレーション
私は COMSOL 希釈種輸送を空間によるシミュレーション
います。先ずは、定義でド
私は COMSOL 希釈種輸送をインポイントプローブを置使います。先ずは、定義でドておきます。リアルタイム
メインポイントプローブを置
拍動の周期的振動を確認でいておきます。リアルタイム
ます(図１７)。ジオメトリ
で拍動の周期的振動を確認では、円のドメインを 4 つにきます(図１７)。ジオメトリ
けます (図１６右下)。初期
では、円のドメインを4 つに件に空間異方性を与えるた

図１７. V の時間変化を計測中に確認。

CO

チュートリアルCOMSOL Multiphysics による計算科学工学 – 生物系（4）

時間分解

が埋めて
ョンが可
の発展に
時間分解

が埋めて
謝辞
の発展に
橋口真宜
図１９. 固体力学(フィジ クス)の線形 弾性材料
で、針による変位を sin(ωt)でテスト。メッシュは
移動メッシュとして、変位に従って動かす。変位
図19
固体力学（フィジクス）
図１９. 固体
力 学 ( フ ィの線形弾性材料で、針による変ジ クス)の線形 弾性材料の絶対量の関数としてIap を増加する。
位を sin(ωt) でテスト。メッシュは移動メッシュとしで、針による変位をsin(ωt)でテスト。メッシュはて、変位に従って動かす。変位の絶対量の関数として
移動メッシュとして、変位に従って動かす。変位
Iap を増加する。
の絶対量の関数として Iap を増加する。

わけます(図１６右下)。初期図17 図１７.
です。また小さい円はペー
V の時間変化を計測中に確認。
V の時間変条件に空間異方性を与えるた化を計測中に確認。
メーカー用です。例えめです。また小さい円はペー
I ap ＝I 0 δ(x, y)のよう

るたびに
謝辞
実験デー
橋口真宜

杉本学君
るたびに
絢子さん
実験デー

新井晋君
杉本学君

絢子さん
参考文献
新井晋君
第１章
[1]
野原
参考文献
2013.
第１章

[1] 野原
第２章～
2013.

スメーカー用です。例えとれば、拍動波面の湧ば出しを限定できます。
I ap ＝I 0 δ(x, y)のよう

図20
左から針で押す。ドメインを弾性体として、変位と
図２ ０.
左か ら針 で押 す。 ドメ イン を弾 性体 とし相関のある Iap を定義すると、刺激に応答して拍動て、変位と相関のある
Iap を定義すると、刺激にが励起、そして伝搬した。

にとれば、拍動波面の湧初期条件に位相差としき出しを限定できます。
図１６右下のようにす

初期条件に位相差としと、半径の長さの波面て図１６右下のようにす
時計回りに回転し、速ると、半径の長さの波面が r に依存しないこと

が時計回りに回転し、速ら、螺旋の波面伝搬がさがr に依存しないことこります(図１８右上)。

図18 V のパルスが波面として伝搬。螺旋（左上）。中心のみ
図１８. V のパルスが波面から、螺旋の波面伝搬がほどの I ap ＝I 0 δ(x, y) Iap ≠0で、中心にペースメーカー（左上）。二つの周として伝搬。螺旋(左上)。
起こります(図１８右上)。
波面の湧き出しを中心 期の異なるドメインは波形の乱れを（下）。
中心のみ
で、中心に
図１８. VIap≠0のパルスが波面
先ほどのI ap ＝I 0 δ(x,y)COMSOL
して円状に波面が広が私は希釈種輸送を使います。先ずは、定義
ペースメーカー(左上)。二として伝搬。螺旋(左上)。
は波面の湧き出しを中心ます。ペースメーカー
でドメインポイントプローブを置いておきます。リア
つ
の周期の
異 なで、
るド
メイ
中心のみ
Iap≠0
中心に
として円状に波面が広がンは波形の乱れを(下)。
ルタイムで拍動の周期的振動を確認できます（図17）。
すね。COMSOL ではδ関
ペースメーカー(左上)。二ります。ペースメーカー
ジオメトリでは、円のドメインを4つにわけます（図16
つの周期の異なるドメイなどの様々な関数に加ンは波形の乱れを(下)。
右下）。初期条件に空間異方性を与えるためです。またですね。COMSOL
ではδ関て、論理式も multiple できるので、さまざまな条
Iap=I0δ(x, y)数などの様々な関数に加を試せます。Iap小さい円はペースメーカー用です。例えば≠0 とする場合は、ドメインを切っのようにとれば、拍動波面の湧き出しを限定できます。
えて、論理式もmultiple
できるので、さまざまな条、式内のパラメータを変えられます。図１８下のよ
初期条件に位相差として図16右下のようにすると、件を試せます。
Iap ≠0 とする場合は、ドメインを切っな、拍動伝搬の速さが異なる任意のドメインも簡単
半径の長さの波面が時計回りに回転し、速さが r に依存て、式内のパラメータを変えられます。図１８下のよ
つなげられて実験条件にあったシミュレーションがしないことから、螺旋の波面伝搬が起こります（図18右
うな、拍動伝搬の速さが異なる任意のドメインも簡単
きます。
上）。先ほどの Iap=I0δ (x, y) は波面の湧き出しを中心とにつなげられて実験条件にあったシミュレーションが
近年では、心臓がメカニカルな動きもフィードバッ
して円状に波面が広がります。ペースメーカーです
できます。
として、拍動を安定化するという予測が話題となっね。COMSOL ではδ関数などの様々な関数に加えて、
近年では、心臓がメカニカルな動きもフィードバッ
います。つまり心臓を弾性体として捉え、伸縮自体論理式も multiple できるので、さまざまな条件を試せまクとして、拍動を安定化するという予測が話題となっ心臓内の電流（電位変化）に寄与するということです。Iap ≠0とする場合は、ドメインを切って、式内のています。つまり心臓を弾性体として捉え、伸縮自体。そこで、COMSOLの versatility と visualization
パラメータを変えられます。図18下のような、拍動伝も心臓内の電流（電位変化）に寄与するということで搬の速さが異なる任意のドメインも簡単につなげられ力を借りて挑戦しました。

す。そこで、COMSOLの versatility と visualizationて実験条件にあったシミュレーションができます。
近年では、心臓がメカニカルな動きもフィードバッ
の力を借りて挑戦しました。
クとして、拍動を安定化するという予測が話題となっています。つまり心臓を弾性体として捉え、伸縮自体
Vol.22, No.4
2017も心臓内の電流（電位変化）に寄与するということです。 そ こ で、COMSOL の versatility と visualization の 力
Vol.22, No.4 2017
を借りて挑戦しました。

計算工学

ョンが可C

応答して拍動が励起、そして伝搬した。
図２ ０. 左か ら針 で押 す。 ドメ イン を弾 性体 とし固体力学のフィジクスで心筋細胞の集合体を線形弾て、変位と相関のある
Iap を定義すると、刺激に固体力学のフィジクスで心筋細胞の集合体を線形弾性性材料としました（図19）。ドメインの動きに連動し応答して拍動が励起、そして伝搬した。
材料としました（図１９）。ドメインの動きに連動して、て、メッシュも移動メッシュの規定メッシュ変位で、メッシュも移動メッシュの規定メッシュ変位で、同じ固体力学のフィジクスで心筋細胞の集合体を線形弾性
同じだけ動かしました。結果が図20です。実験結果のだけ動かしました。結果が図２０です。実験結果の動
材料としました（図１９）。ドメインの動きに連動して、動きに、もっと近づけるためには変位の関数の最適化
きに、もっと近づけるためには変位の関数の最適化がメッシュも移動メッシュの規定メッシュ変位で、同じが必須ですが、実験結果と同様に針刺激で拍動が励起必須ですが、実験結果と同様に針刺激で拍動が励起すだけ動かしました。結果が図２０です。実験結果の動するシミュレーションは出来ました。
るシミュレーションは出来ました。
きに、もっと近づけるためには変位の関数の最適化が4.4 実験結果との考察：細胞の集合体の拍動とメカニ必須ですが、実験結果と同様に針刺激で拍動が励起す
カルな刺激
４．
４ 実験結果との考察：細胞の集合体の拍動とメカるシミュレーションは出来ました。
この心筋細胞の拍動実験と、そのシミュレーション
ニカルな刺激は、それぞれが未完成で、改良の余地があります。しこの心筋細胞の拍動実験と、そのシミュレーション
４．４ 実験結果との考察：細胞の集合体の拍動とメカ
かし、動きのフィードバックのシミュレーションがは、それぞれが未完成で、改良の余地があります。しニカルな刺激COMSOL によって出来そうです。今後は実験系の刺激
かし、動きのフィードバックのシミュレーションがこの心筋細胞の拍動実験と、そのシミュレーション
パラメータの幅を広げて、細胞の応答を観測して、シCOMSOLによって出来そうです。今後は実験系の刺激パ
は、それぞれが未完成で、改良の余地があります。しミュレーションの精度も高めたいです。そして、パララメータの幅を広げて、細胞の応答を観測して、シミ
かし、動きのフィードバックのシミュレーションがメータの生物学的意味としてのチャネルの機能の解明
ュレーションの精度も高めたいです。そして、パラメ
COMSOL によって出来そうです。今後は実験系の刺激パを目指します。
ータの生物学的意味としてのチャネルの機能の解明を
ラメータの幅を広げて、細胞の応答を観測して、シミ目指します。
ュレーションの精度も高めたいです。そして、パラメ
第2章〜第4章のまとめをします。

1.
Enginee
第２章～
Microflu1.
Fluid M
Engine2.
Microflmicroflu

Fluid M
Microflu

2.
423-443
microfl3.
MicroflMicroflu
423-443
scale. R
3.
4.
MicroflElectrooscale. R
microflu
4.
Nanoflu
Electro5.
microflMarzialNanofluT.; Di V
5.
X.; Jova
MarziaLing, X
T.; Di V
Oliver,X.; Jova
Riehn,
Ling, X
WanunuOliver,potentia
Riehn,
sequenc
Wanun
1146-53
ータの生物学的意味としてのチャネルの機能の解明を実験的に考察しやすいようにシンプルな系をつくっ
potentiT., Dir
第２章～第４章のまとめをします。
目指します。
て実験をしているつもりですが、形状、拍動間隔など
sequenmolecul
実験的に考察しやすいようにシンプルな系をつくって
1146-53
実 験 的 に 制 御 す る の は 難 し い で す。 し か し な が ら
nano 20
T., Dir
COMSOL により計測がユーザーフレンドリーになりま
実験をしているつもりですが、形状、拍動間隔など実第２章～第４章のまとめをします。
Kato, Y.
molecuした。COMSOLの特徴である豊富なフィジクスの追加
験的に制御するのは難しいです。
しかしながらCOMSOL
実験的に考察しやすいようにシンプルな系をつくってDNA mnano 2とメッシュの最適化、そして、実験条件の変化に応じにより計測がユーザーフレンドリーになりました。
near an
実験をしているつもりですが、形状、拍動間隔など実Kato, Yた柔軟なパラメータの変化など、実験に忠実なシミュ
controll
COMSOLの特徴である豊富なフィジクスの追加とメッ験的に制御するのは難しいです。しかしながら COMSOL
DNA m
レーションが可能になりました。実験の場合、観測に065502.
シュの最適化、そして、実験条件の変化に応じた柔軟により計測がユーザーフレンドリーになりました。
near an
空間・時間分解能の制限がかかりますが、その隙間も6.
なパラメータの変化など、実験に忠実なシミュレーシ
control
COMSOL の特徴である豊富なフィジクスの追加とメッNanopoCOMSOL が埋めてくれます。COMSOL は、私にとって065502シュの最適化、そして、実験条件の変化に応じた柔軟
今後の研究の発展に欠かせない存在になりました。
6.
なパラメータの変化など、実験に忠実なシミュレーシ
計算工学(??-10)Nanopo

計算工学（22-9）

Vol.23, No.1 2018(??-10)

チュートリアルCOMSOL Multiphysics による計算科学工学 – 生物系（4）謝辞
橋口真宜氏には、新たな系のシミュレーションを始めるたびにアドバイスをいただいていております。また、実験データを献身的に出してくれた石田研太郎博
士、杉本学君、加藤佑太君、坂下直人君、港聖也君、
山下絢子さん、田村圭さん、露木彩葉さん、上原貴宏
君、新井晋君に感謝します。
■参考文献
[1]
[2]

野原勉、応用微分方程式講義、東京大学出版会、2013.
Stone, H. A.; Stroock, A. D.; Ajdari, A., Engineering ﬂows insmall devices: Microﬂuidics toward a lab-on-a-chip. Annu Rev
Fluid Mech 2004, 36, 381-411.

[3]

Sun, T.; Morgan, H., Single-cell microfluidic impedance
cytometry: a review. Microfluidics and Nanofluidics 2010, 8(4), 423-443.
Squires, T. M.; Quake, S. R., Microﬂuidics: Fluid physics atthe nanoliter scale. Rev Mod Phys 2005, 77 (3), 977-1026.
Wang, X.; Cheng, C.; Wang, S.; Liu, S., Electroosmotic pumps
and their applications in microfluidic systems. Microfluidics

[4]
[5]

[6]

[7]

[8]

and Nanofluidics 2009, 6 (2), 145.
Branton, D.; Deamer, D. W.; Marziali, A.; Bayley, H.; Benner,
S. A.; Butler, T.; Di Ventra, M.; Garaj, S.; Hibbs, A.; Huang,X.; Jovanovich, S. B.; Krstic, P. S.; Lindsay, S.; Ling, X. S.;
Mastrangelo, C. H.; Meller, A.; Oliver, J. S.; Pershin, Y. V.;
Ramsey, J. M.; Riehn, R.; Soni, G. V.; Tabard-Cossa, V.;
Wanunu, M.; Wiggin, M.; Schloss, J. A., The potential andchallenges of nanopore sequencing. Nature biotechnology
2008, 26 (10), 1146-53;
Ando, G.; Hyun, C.; Li, J.; Mitsui, T., Directly observing the
motion of DNA molecules near solid-state nanopores. ACS
nano 2012, 6 (11), 10090-7;
Sugimoto, M.; Kato, Y.; Ishida, K.; Hyun, C.; Li, J.; Mitsui, T.,
DNA motion induced by electrokinetic ﬂow near an Au coatednanopore surface as voltage controlled gate. Nanotechnology
2015, 26 (6), 065502.

Google Scholar で検索すると、Nanopore AND FEM(1220)< Nanopore AND COMSOL(1670) と件数が出ます。
[10] Wong, C. T.; Muthukumar, M., Polymer capture by electroosmotic ﬂow of oppositely charged nanopores. J Chem Phys
2007, 126 (16), 164903.
[9][11] http://www.kesco.co.jp/comsol/faq/comsol_geometry_mesh_
post_v51.pdf を参考に
[12] Laohakunakorn, N.; Ghosal, S.; Otto, O.; Misiunas, K.;
Keyser, U. F., DNA interactions in crowded nanopores. Nano
Lett 2013, 13 (6), 2798-802.

計算工学[13] Mayerson, P. L.; Fallon, J. F., The spatial pattern and temporalsequence in which feather germs arise in the white Leghorn
chick embryo. Dev Biol 1985, 109 (2), 259-67;[14] Painter, K. J.; Hunt, G. S.; Wells, K. L.; Johansson, J. A.;
Headon, D. J., Towards an integrated experimental-theoretical
approach for assessing the mechanistic basis of hair andfeather morphogenesis. Interface focus 2012, 2 (4), 433-50.
[15] Ishida, K.; Mitsui, T., Generation of bioengineered feather
buds on a reconstructed chick skin from dissociated epithelial
and mesenchymal cells. Development, growth & differentiation
2016, 58 (3), 303-14.
[16] Just Submitted.
[17] Turing, A. M., The Chemical Basis of Morphogenesis.
Philosophical Transactions of the Royal Society of London.
Series B, Biological Sciences 1952, 237 (641), 37.
[18] Meinhardt, H., Biological pattern formation: new observations
provide support for theoretical predictions. BioEssays : news
and reviews in molecular, cellular and developmental biology
1994, 16 (9), 627-32.
[19] Kondo, S.; Miura, T., Reaction-diffusion model as a framework
for understanding biological pattern formation. Science 2010,
329 (5999), 1616-20.
[20] Maini, P. K.; Woolley, T. E.; Baker, R. E.; Gaffney, E. A.; Lee,S. S., Turing's model for biological pattern formation and therobustness problem. Interface focus 2012, 2 (4), 487-96.
[21] Michon, F.; Forest, L.; Collomb, E.; Demongeot, J.; Dhouailly,D., BMP2 and BMP7 play antagonistic roles in feather
induction. Development 2008, 135 (16), 2797-805.
[22] Hodgkin, A. L.; Huxley, A. F., A quantitative description of
membrane current and its application to conduction andexcitation in nerve. The Journal of Physiology 1952, 117 (4),500-544.
[23] Luo, C. H.; Rudy, Y., A dynamic model of the cardiac
ventricular action potential. I. Simulations of ionic currents
and concentration changes. Circulation Research 1994, 74 (6),1071.
[24] Luo, C. H.; Rudy, Y., A dynamic model of the cardiac
ventricular action potential. I. Simulations of ionic currents
and concentration changes. Circulation research 1994, 74 (6),1071-96.
[25] Fitzhugh, R., Impulses and Physiological States in Theoretical
Models of Nerve Membrane. Biophysical journal 1961, 1 (6),445-66.

（22-10）Vol.23, No.1 2018