さて出だしはカッコいいのですが、かなりトホホなお話をいたします（笑）。故あってある種の偏微分方程式を解いてそれをもとにシミュレーションをやることになりました。このタイプの偏微分方程式は一般解がすでに知られていますので、初期条件と境界条件を入れていくつかの定数を求めます。ところが組んだモデルはそうややこしいものではないのですが、６個の定数がいかさまややこしい式で表されます。決して難しいわけではないのですがあまりに煩雑で計算があっているか間違ってるかがわかりません。実際に数値例をいれてエクセルで計算させた後、おかしいかどうかがわかります。あまりに現実と違う値になったり、グラフを描かせてみるととんでもないことになっているから・・・（笑）。
というわけで初めて使ってみたのが数式処理システムのmaximaというフリーウェアです。６元の連立方程式を入力すると文字式の形で解いてくれます。定積分も積分区間が文字で表されてもお構いなしに数学的に正しい文字式で解を出してきます。恐ろしく強力なシステムで感動すら覚えます（笑）。実はそれまで手計算でシコシコ何カ月もやっていたというあまり笑えないお話です。因みに日本語対応もしているようなので興味のある方はぜひ・・・大変素晴らしいシステムだと思います。
2層膜（気相と液相での拡散）：ラプラス次元ｓでの解析解の導出
*特に断らない限り単位面積当たりとする。
Dsc 角質層拡散係数
Lsc 角質層膜厚さ
Ksc-g 気相-角質層膜間の分配係数
Dve viable skin拡散係数
Lve viable skin厚さ
Kve-g 気相-viable skinの分配係数
hm: 対流伝達係数（気相側）
Cla: 気相中濃度Cg(t)のラプラス変換（一定の時はC0/s）
Acc(s):累積透過量のラプラス変換
x: 各層の皮膚表面方向からの距離
s: ラプラス次元変数（複素数）
解析解記述簡略化のための定数
ds=sqrt(s/Dsc)*Lsc (Dsc=s/(ds/Lsc)^2) (sqrt:square root)
dd=sqrt(s/Dve)*Lve (Dve=s/(dd/Lve)^2)
ラプラス次元での拡散方程式の一般解
角質層(-Lsc =<x=<0)：Csc(s,x)=C*cosh(ds/Lsc*x)+D*sinh(ds/Lsc*x)
Viable skin (Lve =>x=>0)：Cve(s,x)=A*cosh(dd/Lve*x)+B*sinh(dd/Lve*x)
A,B,C,Dは定数
支配方程式(4-7式)の両辺をラプラス変換して
C*s/(ds/Lsc)*sinh(-ds)+D*s/(ds/Lsc)*cosh(-ds)=hm/Ksc-g*(C*cosh(-ds)+D*sinh(-ds))- hm*Cla・・・・（1：気相と角質の境界,7式）
注）符号に注意
A*cosh(dd)+B*sinh(dd)=0・・・（2：シンク条件,4式）
A=Kve-g/Ksc-g*C・・・（3：角質とviable skinとの分配係数,5式）
D*s/(ds/Lsc)=B*s/(dd/Lve)・（4：角質とviable skinの間のflux,6式）
Acc=-(A/(dd/Lve)*sinh(dd)+B/(dd/Lve)*cosh(dd))
（5：Acc=- Dve*∂Cve(s,x)/∂x /s (x=Lve:積分則)）
1）から4）の連立方程式を解いてA,B,C,Dを求めAccを求める。
cosh & sinh：双曲線関数として
Acc(s)=hm*Cla/(gs+gh)
gs=(Ksc-g*Lsc*dd/(Kve-g*Lve*ds)*sinh(dd)*sinh(ds)+cosh(dd)*cosh(ds))*s
gh=(ds/(Ksc-g*Lsc)*cosh(dd)*sinh(ds)+dd/(Kve-g*Lve)*sinh(dd)*cosh(ds))*hm
文献2の拡散方程式(再掲)
1) ∂Csc/∂t ＝Dsc＊∂2Csc/∂x2 (Lve＜ｘ＜Lsc+Lve)
2) ∂Cve/∂t ＝Dve＊∂2Cve/∂x2 (0＜ｘ＜Lve)
3) Cve＝Cv0(x), Csc=Cs0(x) (0<x<Lsc+Lve ,t=0)
4) Cve=0 (x=0, t>0)
5) Csc/Ksc_g=Cve/Kve_g (x=Lve, t>0)
6) Dsc＊∂Csc/∂x= Dve＊∂Cve/∂x (x=Lve,t>0)
7) Dsc＊∂Csc/∂x= hm＊(Cg(t)－Csc/Ksc_g) ( x=Lsc+Lve,t>0)
若干の訂正
x=0を角質層とviable skinの境界にする。
正負の方向を皮膚表面から毛細血管にする。
基本的には3式は必要がない。ラプラス変換では時間無限大まで積分するので・・
*計算上の対応
暴露濃度C0で暴露し続ける時の累積透過量をAcc(t):(t>0)とすると、a時間後からb時間だけC0で暴露した場合の累積透過量は次のようにあらわされる。
Cla=C0*exp(-a*s)/s*(1-exp(-b*s))
Acc(t)=0 :(a>=t)
Acc(t-a) :(a+b>t>a)
Acc(t-a)-Acc(t-a-b):(t>=a+b)
C0に関する比例計算と時間軸のスライドのみで様々なスケジュールのステップ関数の暴露パターンに対応が可能。
https://www.jeea.or.jp/course/contents/01138/
*2層膜の定常状態でのラグタイムLTの算出法
x=Dsc/Lsc^2
y=Ksc*Lsc
z=Dve/Lve^2
u=Kve*Lve
LT=(1/x+1/z)/2-(1/3)*1/(1/(x*y)+1/(z*u))*(y/(x*y)^2+u/(z*u)^2)
------------------------------------------------------------------------------------
まとめ：ラプラス次元での解析解と定数
Acc(s)=hm*Cla/(gs+gh)
gs=(Ksc*Lsc*dd/(Kve*Lve*ds)*sinh(dd)*sinh(ds)+cosh(dd)*cosh(ds))*s
gh=(ds/(Ksc*Lsc)*cosh(dd)*sinh(ds)+dd/(Kve*Lve)*sinh(dd)*cosh(ds))*hm
A=(Cla*Ksc*Kve*Lsc*dd*sinh(dd)*ds*hm)/us
B=-(Cla*Ksc*Kve*Lsc*dd*cosh(dd)*ds*hm)/us
C=(Cla*Ksc*Ksc*Lsc*dd*sinh(dd)*ds*hm)/us
D=-(Cla*Ksc*Kve*Lve*cosh(dd)*ds^2*hm)/us
Csc=-(hm*(Cla*Ksc*Kve*Lve*cosh(dd)*ds^2*sinh((ds*x)/Lsc)-Cla*Ksc*Ksc*Lsc*dd*sinh(dd)*ds*cosh((ds*x)/Lsc)))/us
Cve=(Cla*Ksc*Lsc*ds*hm*(Kve*dd*sinh(dd)*cosh((dd*x)/Lve)-Kve*dd*cosh(dd)*sinh((dd*x)/Lve)))/us
us=((Ksc*Ksc*Lsc^2*dd*sinh(dd)*sinh(ds)+Ksc*Kve*Lsc*Lve*cosh(dd)*ds*cosh(ds))*s+
(Kve*Lve*cosh(dd)*ds^2*sinh(ds)+Ksc*Lsc*dd*sinh(dd)*ds*cosh(ds))*hm)
気相と角質層境界での累積透過量xst
xst=(Cv/s**Ksc*Lsc^2*dd*sinh(dd)*hm)*sinh(ds)/us+(Cv/s*Ksc*Lsc*Kve*Lve*cosh(dd)*ds*hm)*cosh(ds)/us
角質層とviable skin境界での累積透過量vexst
vexst=(Cv/s*Ksc*Lsc*Kve*Lve*cosh(dd)*ds*hm)/us
角質層内薬物量Msc
Msc=xst-vexst
viable skin内薬物量Mve
Mve=vexst-Acc
