死後経過時間推定するために方程式をscipy.optimizeを利用して解く

どうして死後経過時間を求めるのか

(引用;名探偵コナン「甘く冷たい宅配便」)
コナン君が遺体を見つけた時、かなり最初の段階で死後経過時間について考えますよね？死亡推定時刻が分かれば大きく容疑者を絞り込むことができますしね＼(^o^)／
そこで今回は死後経過時間について詳しい求め方を紹介します！

例題

室温16℃に設定された冷蔵車に男性の遺体があった。男性の体重は86kg、コナン君が遺体の直腸内温度を計ったところ、27℃であった。男性の死後経過時間を求めよ。

ニュートンの冷却の法則

まず、簡略化した事象を考えてみましょう。コップに入れられた熱いコーヒーが周囲の温度（空気）によって冷まされる問題を考えます。この時、次の「ニュートンの冷却の法則」が近似的に成り立つと言われています。
ニュートンの冷却の法則を表す微分方程式は

$\frac{dT}{dt} = -γ(T - T_s)$ (1)

で表されます。ここで、各文字は
$T$：湯の温度
$T_s$：周囲の温度
$\gamma$：冷却定数
を表しています。この微分法方程式は解析的に解くことができます。まず、(1)式の両辺を$T-T_s$で割ります。
$\frac{1}{T-T_s}dT = -γdt$

両辺を積分して積分定数をCとすると、
$log (T - T_s) = -γt + C$
となります。この式をTについて解くと下式のようになります。
$T = T_s + e^{-γt + C}$ (2)
初期条件$T(0)=T_0$を課すと、積分定数Cは
$C = log (T_0-T_s)$

となります。これを(2)式に代入すれば、解析解は以下で与えられます。
$T(t) = T_s + (T_0 - T_s) e^{-γt}$

（また、これを変形して）
$\frac{T(t)-T_s}{T_0-T_s} = e^{-γt}$

直腸温の連続測定による死後経過時間の推算

次に、ニュートンの冷却の法則を利用して死後経過時間を求めます。「直腸温の連続測定による死後経過時間の推算（早大院先進理工）」によると遺体の場合、微分方程式は次の様に変形されます。

\frac{T_r-T_e}{T_o-T_e}=1.25e^{B*t}-0.25e^{5B*t}\qquad(T_e\le23.2)\\

\frac{T_r-T_e}{T_o-T_e}=1.11e^{B*t}-0.11e^{10B*t}\qquad(T_e\ge23.2)

ただし、
$T_r$ は直腸温 ℃
$T_e$ は環境温 ℃
$T_o$ は死亡時の直腸温であり 37.2℃ とする
$T_e$ は環境温 ℃ で測定時間の平均値で一定とする
$W$ 死体の体重kg
$b=1$ は補正係数
$B$は Newtonの冷却定数
$t$ は死後経過時間(h)
とします。
これを元に早速今回の例題を解きましょう！

.py

import sympy

def death_time(Tr,Te,W):
    # Tr: 直腸温,Te, 環境温, W: 体重
    To = 37.2 # 死亡時直腸温
    b = 1 # 補正係数
    E = sympy.S.Exp1
    B = sympy.Symbol('B') # 変数B: Newtonの冷却定数
    t = sympy.Symbol('t') # 変数t: 死後経過時間
    
    # 左辺が0になるように方程式を連立する
    equation1 = 1.25*(E**(B*t)) - 0.25*(E**(5*B*t)) - (Tr-Te)/(To-Te)
    equation2 = 1.11*(E**(B*t)) - 0.11*(E**(10*B*t)) - (Tr-Te)/(To-Te)
    equation3 = -1.2815*((b*W)**(-0.625)) + 0.0284 - B # 冷却定数の式

    if Te <= 23.2:
        # 環境温が23.2℃以下の場合
        return sympy.solve([equation1, equation3])
    else:
        # 環境温が23.2℃より高い場合
        return sympy.solve([equation2,equation3])

print(death_time(Tr=27, Te=16, W=86)) # 今回の問題の設定       

$ $参考文献によると$T_e$は環境温の平均を使いますが、今回は発見時の室温の16℃を$T_e$としました。
結果は言わなくても分かると思いますが、scipyでこんな累乗変数を含む計算を直接できるわけもなく、朝起きても計算を続けていました(泣)。参考文献通りに試行錯誤法の考えの元で解くことにしました。

(引用Krylov Subspace Accelerated Newton Algorithm: Application to Dynamic Progressive Collapse Simulation of Frames)

関数の微分を求め方程式を反復の都度解く処理が大きい為、上図の様に微分を差分で近似させたクリロフ部分交換法を用いて計算することにしました。

.py

import math
import sympy
from scipy.optimize import newton_krylov

#(3) 式に相当する。体重 (kg) からNewtonの冷却定数Bを算出。To (死亡時直腸温) は37.2とする
def predict_B(W):
    B = sympy.Symbol('B') # 変数: ニュートンの冷却定数
    b = 1 # 補正係数
    equation3 = -1.2815*((b*W)**(-0.625)) + 0.0284 - B
    return sympy.solve([equation3])[B]

# (1), (2) 式に相当。環境温 Te で場合分け
def predict_deathtime(W, Tr, Te, ta):
    # 求めるものはt. オプティマイザを利用する上での t の初期値を ta とする.
    Bx = predict_B(W) # 求値したB
    E = sympy.S.Exp1
    To = 37.2 # 死亡時直腸温
    
    def F(x):
        # F(x) = 0 の解を求めることができる
        if Te <= 23.2:
            return 1.25*(E**(Bx*x)) - 0.25*(E**(5*Bx*x)) - (Tr-Te)/(To-Te)              
        else:
            return 1.11*(E**(Bx*x)) - 0.11*(E**(10*Bx*x)) - (Tr-Te)/(To-Te)

    sol = newton_krylov(F, ta, method='lgmres')
    return sol

# 死後経過時間の初期値を0にすると微分値が求められないので注意する

predict_deathtime(86, 27, 16, 2)

17.19191007875492

よって、今回の例題では17時間ほど前に死亡したことが分かりました！

参考文献

計算科学をはじめよう！ ニュートンの冷却法則①
直腸温の連続測定による死後経過時間の推算（早大院先進理工）○（学）井上 幹康*（正）酒井 清孝（早大高等研）（正）山本 健一郎（防衛医大）（正）金武 潤