Okay, so with the approval of @charlesm93 Iâ€™ll post the ODE describing the pharmacokinetic model. Iâ€™ll tack on a short derivation from first principles for good measure and Iâ€™ll tag people that probably could tell me at once if it is absurd: @wds15, @yizhang, @andrewgelman

Itâ€™s actually not *that* wild.

First: In addition to the parameters in Table 1 of the original paper, there are four additional parameters that are measured and assumed to be known â€śexactlyâ€ť:

- The lean body mass [kg],
- The mass fraction of fat [N/A],
- The pulmonary volume flow [l/min],
- The concentration of PERC in the air during exposure [PPM] with 1PPM=409/144 mug/l.

From the lean body mass [kg] and the mass fraction of fat [N/A] we can compute the total volume of fat [l] from the density of fat (.92 g/ml). The body of each person gets modeled as four distinct compartments which are connected with each other via the circulatory system, mixing the blood â€ścoming outâ€ť of each compartment, then passing through the lungs and mixing with in/exhaled air, and then again feeding the blood into each compartment. As described in the original paper, section 2.1, the concentration C_s [mg/l] in each compartment follows the ODE

dC_s/dt=(C_{art}-C_s/P_s)F_s/V_s,

except for the liver compartment where we have the additional mass sink

V_{max}C_s/(V_s(K_m+C_s)).

The latter term is just the Michaelisâ€“Menten ansatz, while the first term derives from the conservation of mass for each compartment. Namely, for the mass in each compartment M_s=C_sV_s [mg] from mass conservation (rate of change in mass = mass flow in - mass flow out) we have

dM_s/dt=(C_{art}-C_{out,s})F_s

where F_s [l/min] is the blood volume flow per compartment, C_{art} [mg/l] is the PERC concentration in the artertial blood (coming from the lungs) and C_{out,s} [mg/l] is the PERC concentration in the blood exiting the compartment. Assuming equilibrium at the compartment exit we get

C_{out,s}=C_s/P_s,

which together with the mass balance yields the original ODE.

There are now still three parts missing. We still need

- an expression for C_{art} [mg/l] (the PERC concentration in the blood exiting the lungs),
- an expression for C_{venous} [mg/l] (the PERC concentration in the mixed blood leaving the compartments and entering the lungs) and
- an expression for C_{exhaled} [mg/l] (the PERC concentration in the exhaled air).

The latter two quantities get measured (with lognormal noise).

The simplest part is the expression for C_{venous}. This is just the concentration of the mixed blood flows exiting the compartments, i.e.

C_{venous} = \sum_s C_{out,s}F_s/F_{venous}

with the total blood flow F_{venous} = \sum_s F_s.

The concentration in the exhaled air can be derived from the assumption of equilibrium between â€śalveolar air [â€¦] and arterial bloodâ€ť (section 2.1). Therefore we just get

C_{exhaled} = C_{art} / P_{blood/air}.

Finally, C_{art} we get again from conservation of mass. The mass flow of PERC entering the lungs has to be equal to the mass flow exiting the lungs. Crucially, there are two ways that PERC can enter and exit the lungs, via the inhaled air+the incoming venous blood or via the exhaled air+the outgoing venous blood. The mass flow [mg/min] is just the concentration [mg/l] times the volume flow[l/min], hence we get

C_{inhaled}F_{alveolar}+C_{venous}F_{venous}=C_{exhaled}F_{alveolar}+C_{art}F_{art}.

With F_{venous}=F_{art} (volume flow of PERC is negligible compared to volume flow of blood) and plugging in C_{exhaled} = C_{art} / P_{blood/air} we get

C_{art} = \frac{C_{inhaled}F_{alveolar}+C_{venous}F_{venous}}{F_{alveolar}/P_{blood/air}+F_{venous}}.

There are actually still some things missing.

First, F_{alveolar} gets computed from the given pulmonary volume flow via the magical constant .7 as F_{alveolar}=.7F_{pulmonary}.

Furthermore, F_{venous} gets computed using the Ventilation/perfusion ratio (VPR) from F_{alveolar}

as F_{venous}=F_{alveolar}/(VPR).

Stay tuned or not while I update this post.

There are two more magical things we have to know before we can hopefully reproduce the results from the paper.

- Apparently, the V_{max} [mg/min] from the equations is not the maximal metabolic rate from table 1. Instead,
`VMI`

from table 1 has units [mg/min / kg^(-.7)], such that

V_{max} = \mathrm{VMI}*\mathrm{bodyWeight}^{0.7}.
- The relative volumes of the four compartments excepting the fat sum to a magically precise .837 (relative to the lean body volume).

The first point makes me wonder whether Iâ€™ve got the units right for `VMI`

([mg/min]?) and `KMI`

([mg/l] to match the `C_s`

?).

As for reproducing the results from the paper, weâ€™ll not do this (yet). For now, consider that

and that

Before we try to fit the full model it might be helpful to look at the ODE a little bit more closely. The whole ODE reads

dC_s/dt=(C_{art}-C_s/P_s)F_s/V_s-\delta_{s,4} V_{max}C_s/(V_s(K_m+C_s))

where s runs from 1 to 4 and \delta_{i,j}=\mathrm{int}(i==j).

There are two things which are not immediately obvious from this formulation. Neglecting the Michaelis-Menten term

- this is a
*linear* ODE and
- there should be a perfect anticorrelation between the volumes V_s and the partition coefficients P_s.

**Claim 1:**

The only term in the ODE without the Michaelis-Menten term that could potentially be nonlinear is C_{art}. However, itâ€™s not:

C_{art}=\frac{C_{inhaled}F_{alveolar}+C_{venous}F_{venous}}{F_{alveolar}/P_{blood/air}+F_{venous}}

The only term that depends on our time-dependent state C_s is C_{venous}, which in turn is just

C_{venous}=\sum_s C_{out,s}F_s/F_{venous}=\sum_k C_{k}F_k/(F_{venous}P_k).

All in all we get

dC_s/dt = (\frac{C_{inhaled}F_{alveolar}+\sum_k C_{k}F_k/P_k}{F_{alveolar}/P_{blood/air}+F_{venous}}-C_s/P_s)F_s/V_s

which is linear in the state C_s. We may write the ODE as

dC_s/dt = \sum_k A_{sk}C_k+f_s

with constant source term (which is zero after exposure)

f_s = \frac{C_{inhaled}F_{alveolar}}{F_{alveolar}/P_{blood/air}+F_{venous}}F_s/V_s

and matrix

A_{sk}=\frac{F_s}{V_s}(\frac{F_k}{P_k(F_{alveolar}/P_{blood/air}+F_{venous})}-\frac{\delta_{s,k}}{P_s}).

**Claim 2**:

The second claim becomes clear once we reformulate the ODE in terms of C_{out,s}=C_s/P_s. We get

dC_{out,s}/dt=dC_s/dt/P_s = \sum_k B_{sk}C_{out,k}+g_s

with new source term

g_s = f_s/P_s = \frac{F_s}{V_s P_s}\frac{C_{inhaled}F_{alveolar}}{F_{alveolar}/P_{blood/air}+F_{venous}}

and new matrix

B_{sk}=A_{sk}\frac{P_k}{P_s}=\frac{F_s}{V_s P_s}(\frac{F_k}{F_{alveolar}/P_{blood/air}+F_{venous}}-\delta_{s,k}).

Stay tuned for the next part where we reveal

- the exact solution of the linear ODE
- an attempt at fitting the linear ODE with the original parametrization
- a hopefully successful attempt at fitting a reduced linear ODE
- a tool to quickly explore the full, nonlinear ODE and
- a hopefully successful and quick fit of a reduced nonlinear ODE.