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.