Unsupervised learning of pharmacokinetic responses

Abstract

Pharmacokinetics (PK) is a branch of pharmacology dedicated to the study of the time course of drug concentrations, from absorption to excretion from the body. PK dynamic models are often based on homogeneous, multi-compartment assumptions, which allow to identify the PK parameters and further predict the time evolution of drug concentration for a given subject. One key characteristic of these time series is their high variability among patients, which may hamper their correct stratification. In the present work, we address this variability by estimating the PK parameters and simultaneously clustering the corresponding subjects using the time series. We propose an expectation maximization algorithm that clusters subjects based on their PK drug responses, in an unsupervised way, collapsing clusters that are closer than a given threshold. Experimental results show that the proposed algorithm converges fast and leads to meaningful results in synthetic and real scenarios.

Appendices

Appendix 1: Experimental setup

As described in Sect. 3.1, the M-step maximizes the likelihood over the parameters $\beta _{1\ell }$ and $\beta _{2\ell }$, for each cluster $\ell$, resorting to a numerical method; indeed, it is overwhelming to obtain the analytical solution of the transcendental system of equations that provided the maxima. Next, we detail a simple optimization and stopping criteria for coordinate descent method used in the experiments presented in Sect. 4.

1.1 Optimization

To improve the convergence rate of coordinate descent method we made an improvement in the method by intercalating the (one-dimensional) Newtown iterations for variables $\beta _{1\ell }$ and $\beta _{2\ell }$, thus, solving Eqs. (4) and (5) simultaneously. With such modification, the method converged significantly faster, without deteriorating the results. Recall that Newton(b, h) is an iterative method to find a root of a function h given an initial approximation b, such that, in the d-th iteration of Newton’s method we have

$$\begin{aligned} b^{(0)}=b\text { and }b^{(d+1)}=b^{(d)}-\frac{h(b^{(d)})}{h'(b^{(d)})}. \end{aligned}$$

We can recast the M-step for maximizing $\beta _{1\ell }$ and $\beta _{2\ell }$ as the problem of finding $b_1$ that maximizes $h_1(b_1,b_2)$, given $b_2$ fixed, and $b_2$ that maximizes $h_2(b_1,b_2)$, given $b_1$ fixed, simultaneously. A simple approach to address this problem is to interleave the iterations of Newton’s method as $b_1^{(0)}, b_2^{(0)}, b_1^{(1)}, b_2^{(1)},\dots$ so that

$$\begin{aligned} b_1^{(d+1)}=b_1^{(d)}-\frac{h_{1}(b_{1}^{(d)},b_{2}^{(d)})}{\frac{\partial h_{1}}{\partial b_1}(b_{1}^{(d)},b_{2}^{(d)})} \text { and } b_2^{(d+1)}=b_2^{(d)}-\frac{h_{2}(b_{1}^{(d)},b_{2}^{(d)})}{\frac{\partial h_{2}}{\partial b_2}(b_{1}^{(d)},b_{2}^{(d)})}. \end{aligned}$$

We took this approach where

$$\begin{aligned} b_1^{(0)}=\beta ^{(k)}_{1\ell }\text { and }b_2^{(0)}=\beta ^{(k)}_{2\ell }. \end{aligned}$$

Note that k is used for the k-th iteration of the EM algorithm, whereas d (in the above case $d=0$) is used for the d-th iteration of the Newton’s method. So the previous estimates of $\beta _{1\ell }$ and $\beta _{2\ell }$ by the EM algorithm, are the initial approximation for the roots in Newton’s method. In this case, we let

$$\begin{aligned} h_1(b_1,b_2)=h^{(k)}_{1\ell }(b_1,b_2)\text { and }h_2(b_1,b_2)=h^{(k)}_{2\ell }(b_1,b_2), \end{aligned}$$

with $h^{(k)}_{1\ell }(b_1,b_2)$ defined as

$$\begin{aligned} \displaystyle \sum _{i=1}^N\sum _{j=1}^n\Big (-\frac{\alpha ^{(k+1)}_\ell }{v^{(k)}_\ell }X_{i\ell }^{(k)}t_je^{-b_1t_j} \left( y_{ij}-\alpha _\ell ^{(k+1)}(e^{-b_1t_j}-e^{-b_2t_j})\right) \Big ), \end{aligned}$$

and $h^{(k)}_{2\ell }(b_1,b_2)$ defined as

$$\begin{aligned} \displaystyle \sum _{i=1}^N\sum _{j=1}^n\Big (\frac{\alpha _\ell ^{(k+1)}}{v^{(k)}_\ell }X_{i\ell }^{(k)}t_je^{-b_2t_j} \left( y_{ij}-\alpha _\ell ^{(k+1)}(e^{-b_1t_j}-e^{-b_2t_j})\right) \Big ). \end{aligned}$$

Observe that functions $h^{(k)}_{1\ell }$ and $h^{(k)}_{2\ell }$ as defined in Sect. 3.1 have only one argument, while here they have two for the purpose of applying the envisage optimization. With this optimization we start with $b^{(0)}_1=\beta ^{(k)}_{1\ell }$ and $b^{(0)}_2=\beta ^{(k)}_{2\ell }$, then we iterate the Newton’s method until convergence in, say, d iterations, and finally, we update $\beta _{1\ell }$ and $\beta _{2\ell }$ as $\beta ^{(k+1)}_{1\ell }=b^{(d)}_1$ and $\beta ^{(k+1)}_{2\ell }=b^{(d)}_2$.

1.2 Stopping criterion

In our implementation, the $\text {Newton}(b,h)$ method is stopped at iteration d if

$$\begin{aligned} \mid b^{(d)}-b^{(d-1)} \mid <10^{-10}. \end{aligned}$$

Moreover, if the number of iterations exceeds $10^4$ the method is aborted (with a warning to the user). However, in practice, this iteration limit was never achieved in the experiments performed in Sect. 4.

In our particular setting, as we are computing at each step both $\beta _{1\ell }$ and $\beta _{2\ell }$, we stop when both approximations fulfill the previous criteria.

1.3 Knee analysis

For the clustered curves to have biological meaning (e.g., concentrations need to be positive, rates cannot become exponentially faster, etc.), we imposed that

$$\begin{aligned} \beta _{1\ell }, \beta _{2\ell }\in (0,5)\text { and }\beta _{1\ell }\le \beta _{2\ell }. \end{aligned}$$

If Newton’s method ends with $\beta _{1\ell }>\beta _{2\ell }$ then we swap $\beta _{1\ell }$ with $\beta _{2\ell }$; in practice, this never happened in the experiments performed in Sect. 4. Nonetheless, if $\beta _{1\ell }\notin (0,5)$ or $\beta _{2\ell }\notin (0,5)$, we perform a knee analysis to elicit a value within the allowed range.

In a general setting of $\text {Newton}(b,h)$, assume that the method converges in iteration d and $b^{(d)}$ is out of some allowed range. It might be very well the case $b^{(d)}$ is out of range but the images of h within the allowed range are very closed to $h(b^{(d)})$. In this case, the method does not converge inside the range because $b^{(d)}>b^{(d-1)}$ but $h(b^{(d)})\approx h(b^{(d-1)})$. Therefore, the method should be stopped when

$$\begin{aligned} \mid h(b^{(d)})-h(b^{(d-1)}) \mid <10^{-10}, \end{aligned}$$

which is called a knee analysis.

In our particular case, we perform this knee analysis to avoid $\beta _{1\ell } \notin (0,5)$ and $\beta _{2\ell } \notin (0,5)$. If, even performing this analysis, we have $\beta _{1\ell } \notin (0,5)$ then we keep the previous parameters approximation, that is, $\beta _{1\ell }^{(k+1)}=\beta _{1\ell }^{(k)}$, and the EM algorithm proceeds. The same analysis is carried out for $\beta _{2\ell }$.

Appendix 2: Implementation details

In this appendix, we present implementation details related to the proposed EM algorithm.

1.1 Initial number of clusters

The initial number of clusters is an user-defined parameter. However, if this value is not given, we set $M = \frac{N}{3}$.

1.2 Stopping criterion

Usually, an EM algorithm stops when the difference between consecutive values of the parameters reaches some threshold. In our case, the impact of the parameters $\beta _{1\ell }$ and $\beta _{2\ell }$ overwhelms the impact of all remaining parameters in the definition of the clusters; recall that $\beta _{1\ell }$ and $\beta _{2\ell }$ are exponents in Eq. (2). For this reason, only $\beta _{1\ell }$ and $\beta _{2\ell }$ are considered in the stopping criterion of the proposed EM algorithm. Specifically, the EM algorithm stops when

$$\begin{aligned} \frac{\mid \beta _{1\ell }^{(k+1)}-\beta _{1\ell }^{(k)}\mid }{\mid \beta _{1\ell }^{(k+1)}\mid }\le 10^{-6}\text { and } \frac{\mid \beta _{2\ell }^{(k+1)}-\beta _{2\ell }^{(k)}\mid }{\mid \beta _{2\ell }^{(k+1)}\mid }\le 10^{-6}, \end{aligned}$$

for all clusters $\ell \in \{1,\dots ,M\}$.

1.3 Cluster thresholds

The thresholds for disregarding and merging clusters were defined (empirically) as $\overline{W}=0.025$ and $\overline{L}=1$.