Comment on how to computationally solve the bi-exponential formula behind Henssge´s estimation of time since death nomograms [pHp; after ten years]

My Henssge online tool (small layout; large layout) solves the bi-exponential formula as laid out in Henssge and Madea (2004) Estimation of the time since death in the early post-mortem period. Forensic Science International 144 (2004) 167-175.

Henssge at the time still regarded the solution to be obtained, by calculators, to necessarily employ some iteration procedure and I do remember some Casio and Hewlett Packard programs that were programmed with iterations. They took a while to get there.

It is now 2015.

I implemented this pHp online tool roughly ten years ago, basically without thinking a lot, in 2005, because the approach is obvious.

I only sat down to write this because it was stated that this could not be done using pHp, when our IT guys proposed that, so the bet was over 1 cup of coffee, which I won #holdmybeer.

I love pHp because it stays where you put it, and it is superbly accessible. Obviously this works, and it works very fast.

The pHp code was hand written using a simple text editor.

It later was modified to contain value delimiting comments and tables helping with the correction factor and so on. It was also validated, in addition to or own tests it was later validated by South African forensic pathologists.

There is a small layout geared towards mobile device usage and a large layout with larger displays in mind; that also contains ancillary information such as corrective factor tables.

So, how is that done?

In essence, following exactly Henssge, we establish the estimate for the time of death $t$ as a function of body core temperature $T_B$ , environmental temperature $T_U$ , initial body core temperature at time of death $T_D$ (which is usually but not necessarily $37,2^{\circ} C$ ), and as function of the body mass $M$ to be estimated here and the corrective factor $f$ (basically modelling environmental aspects of cooling such as humidity and body covers, see the large layout for detailed numerical suggestions).

So all one wants to do is find the solution for something like this equation:

$\frac{T_B-T_U}{T_D-T_U} = 1.25e^{Zt} - 0.25e^{5Zt}$ (1)

$\frac{T_B-T_U}{T_D-T_U} = 1.11e^{Zt} - 0.11e^{5Zt}$ (2)

with

$Z=(-1.2815 \times (f \times M^{-0.625})+0.0284)$ (3)

Up to an ambient temperature T_U≤23.2^∘C: Eq. 1; for any ambient temperature T_U≥23.3^∘C: Eq. 2.

With pHp, the approach is obviously brute force to be fast, and it seems that I used the thinking behind the Newton Raphson method to pragmatically get the job done.

This equation is not hard to understand: given realistic values for all given variables, you will have estimated or guessed the correct time of death estimate $t$ once the terms on both sides of the equal signs are actually equal, and you are close to guessing the correct time of death estimate according to this model once you approximate equality. So when $X = Y$ (as in Eq. 1 or 2), it follows that $X - Y = 0$ .

Just performing a forward calculation of both terms of that equations 1 and 2, using a large vector of all $N$ possible time of death estimates $t_A(i)$ with $n=\{1,N\}$ will yield $N$ different equations.

The difference between both terms (on the right and left side of the equation) tells you how close you guessed, and you have guessed well once the difference approximates zero, and if we calculate a few hundred options for $X$ and $Y$ across a set of some possible time of death estimates, the whole set of $X$ and $Y$ will contain a best guess i where $X_{i} - Y_{i}$ is minimal, if not zero.

So, I just shove a vector $t_A$ of all possible sensibly relevant times of death $n=\{1,N\}$ (tip: some $N=300$ discrete values are totally enough or so to cover the extent of the method somewhere between time of death estimates ranging from around an hour to two or so days) through the following equation in order to determine, at once, all $N$ values for the difference vector $D$ with $D=\frac{T_B-T_U}{T_D-T_U} -( 1.25e^{(-1.2815 \times (f \times M^{-0.625})+0.0284) \times t_A} - 0.25e^{5 \times (-1.2815 \times (f \times M^{-0.625})+0.0284) \times t_A})$ and as a solution, identify the one value $t_A(i)$ for which $D(i)=min(|D|)$ .

So in essence, within two steps of actually “calculating” something, basically, you are done.

For better method application, the coefficients differ, depending on specifics, maybe see the paper with regard to this [link].