Énoncé du problème

Nous considérons une cascade parent–métabolite pertinente pour la pharmacocinétique/pharmacodynamie (PK/PD) :

\[ A \xrightarrow{k_f} B \xrightarrow{k_e} C, \]

où $A$ est le parent administré (ou prodrogue), $B$ est la moitié active qui entraîne l’effet pharmacodynamique (PD), et $C$ est un produit inactif ensuite éliminé. L’élimination de $C$ n’a pas d’importance pour l’évolution temporelle de $B$ et est ignorée dans toute la suite.

Lien clinique/PD. Supposons que l’effet pharmacodynamique est directement et instantanément entraîné par la concentration du métabolite actif, $E(t)\propto B(t)$. Par conséquent, le temps de l’effet maximal est égal au temps de concentration maximale du métabolite actif, noté $t_\text{peak}$.

Connus

Demi-vie d’élimination du métabolite actif $B$ : $T_{B\to C}>0$, donnant le taux d’élimination du premier ordre $k_e=\ln 2/T_{B\to C}$.
Temps observé jusqu’à la concentration maximale du métabolite actif (et de l’effet) : $t_\text{peak}>0$.
Condition initiale reflétant un bolus IV (dose normalisée) : $A(0)=1,\;B(0)=0$.

L’objectif est de calculer la demi-vie de formation du parent vers le métabolite actif : $T_{A\to B}=\ln 2/k_f$ uniquement à partir des paramètres connus.

Modèle pharmacocinétique

Il s’agit d’une cascade linéaire, à sens unique, du premier ordre sans rétroaction :

Formation de $B$ à partir de $A$ au taux $k_f$ (unités : 1/temps). Pharmacologiquement, $k_f$ capture l’apparition nette de la moitié active à partir de son précurseur (par exemple, activation métabolique), en supposant que la distribution est rapide par rapport à la formation/élimination pour cette analyse simplifiée.
Élimination de $B$ au taux $k_e$ (unités : 1/temps), qui est déterminée par la demi-vie mesurée $T_{B\to C}$.

Les EDO de bilan de masse sont

\[ \dot A(t)=-k_f A(t), \qquad \dot B(t)=k_f A(t)-k_e B(t), \]

avec $A(0)=1$ et $B(0)=0$. Celles-ci impliquent

\[ A(t)=e^{-k_f t}. \]

En substituant dans l’équation de $B$, on obtient une EDO linéaire avec entrée connue $k_f A(t)$. La résolution (via le facteur intégrant) donne, pour $k_f\neq k_e$,

\[ B(t)=\frac{k_f}{k_e-k_f}\Big(e^{-k_f t}-e^{-k_e t}\Big),\qquad t\ge 0. \]

Interprétation. $B(t)$ est la différence de deux exponentielles : elle augmente initialement lorsque la formation domine, puis diminue lorsque l’élimination domine. Le pic se produit lorsque l’entrée et la sortie instantanées de $B$ sont égales.

2) Temps de pic du métabolite actif (et de l’effet)

Différenciez $B(t)$ et posez égal à zéro :

\[ \frac{dB}{dt} =\frac{k_f}{k_e-k_f}\Big(-k_f e^{-k_f t} + k_e e^{-k_e t}\Big)=0 \quad\Longrightarrow\quad k_f e^{-k_f t}=k_e e^{-k_e t}. \]

En prenant les logarithmes naturels :

\[ \ln k_f - k_f t = \ln k_e - k_e t \;\;\Longrightarrow\;\; t_\text{peak}=\frac{\ln(k_f/k_e)}{\,k_f-k_e\,}. \]

Cette formule est dimensionnellement cohérente (temps des deux côtés) et valide pour $k_f\neq k_e$.

Cas spécial $k_f=k_e$. En utilisant la règle de l’Hôpital ou en résolvant l’EDO dégénérée, on trouve

\[ B(t)=k_f t\,e^{-k_f t},\qquad t_\text{peak}=1/k_f=\frac{T_{B\to C}}{\ln 2}, \]

c’est-à-dire que lorsque la formation et l’élimination sont également rapides, le pic se produit à une durée de vie moyenne et $T_{A\to B}=T_{B\to C}$.

3) Inversion de la formule du temps de pic pour obtenir $k_f$

Notre objectif est de calculer $k_f$ (donc $T_{A\to B}$) à partir des $k_e$ connus et du $t_\text{peak}$ observé. Définissons le rapport sans dimension

\[ r \equiv \frac{k_f}{k_e}>0,\qquad a \equiv k_e\,t_\text{peak}=\frac{(\ln 2)\,t_\text{peak}}{T_{B\to C}}>0. \]

La relation du temps de pic devient

\[ a=\frac{\ln r}{r-1}. \]

Cette équation transcendante a une solution analytique utilisant la fonction Lambert $W$. Procédons comme suit :

\[ \begin{aligned} a(r-1)&=\ln r &&\Longrightarrow& e^{a(r-1)}&=r\\ &&&\Longrightarrow& r\,e^{-a r}&=e^{-a}\\ &&&\Longrightarrow& (-a r)\,e^{-a r}&=(-a)\,e^{-a}. \end{aligned} \]

Par définition de la fonction Lambert $W$ ($W(z)\,e^{W(z)}=z$), nous identifions

\[ -a r = W\!\big(-a e^{-a}\big) \quad\Longrightarrow\quad r = -\frac{1}{a}\,W\!\big(-a e^{-a}\big). \]

Récupérez $k_f=r\,k_e$, puis convertissez en la demi-vie souhaitée

\[ \boxed{\;T_{A\to B}=\frac{\ln 2}{k_f} = -\,\frac{(\ln 2)\,t_\text{peak}}{\,W_b\!\big(-a e^{-a}\big)}\;}, \qquad a=\frac{(\ln 2)\,t_\text{peak}}{T_{B\to C}}. \]

Qu’est-ce que Lambert $W$ ? C’est l’inverse multivalué de $x\mapsto x e^{x}$. De nombreuses bibliothèques scientifiques l’implémentent comme lambertw. Dans notre cas, l’argument $z=-a e^{-a}$ se situe dans $(-e^{-1},0)$, où deux branches réelles existent : la branche principale $W_0$ et la branche inférieure $W_{-1}$.

Sélection de branche (pourquoi deux $W$ possibles, et laquelle utiliser ?)

L’application $r\mapsto a=\frac{\ln r}{r-1}$ est monotone décroissante pour $r>0, r\neq 1$, avec les limites

\[ \lim_{r\to 0^+} a=+\infty,\qquad a(1)=1,\qquad \lim_{r\to\infty} a=0^+. \]

Par conséquent :

Si $r<1$ (formation plus lente que l’élimination, $k_f<k_e$), alors $a>1$ (le pic se produit après une durée de vie moyenne de $B$).
Si $r>1$ (formation plus rapide que l’élimination, $k_f>k_e$), alors $a<1$ (le pic se produit avant une durée de vie moyenne de $B$).
Si $r=1$, alors $a=1$ (le cas spécial ci-dessus).

Pour $z=-a e^{-a}\in(-e^{-1},0)$, les deux branches réelles se comportent comme :

$W_0(z)\in(-1,0)$ $\Rightarrow$ donne $r\in(0,1)$.
$W_{-1}(z)\in(-\infty,-1]$ $\Rightarrow$ donne $r\in[1,\infty)$.

Règle (en termes d’observables) :

\[ b=\begin{cases} 0, & a\gt 1 \quad (\text{utiliser } W_0 \Rightarrow k_f \lt k_e \Rightarrow T_{A\to B} \gt T_{B\to C}),\\ -1,& a\lt 1 \quad (\text{utiliser } W_{-1} \Rightarrow k_f \gt k_e \Rightarrow T_{A\to B} \lt T_{B\to C}),\\ \text{l'une ou l'autre},& a=1 \quad (\text{donne } k_f=k_e). \end{cases} \]

Solution analytique finale (prête à l’emploi)

Soit

\[ a=\frac{(\ln 2)\,t_\text{peak}}{T_{B\to C}},\qquad z=-a e^{-a}. \]

Alors

\[ \boxed{\;T_{A\to B} = -\,\frac{(\ln 2)\,t_\text{peak}}{\,W_b(z)\,}\;}, \]

avec la branche $b$ choisie via la règle ci-dessus (comparez $a$ avec 1).

Vérification des unités. $a$ et $z$ sont sans dimension. Le numérateur a des unités de temps, $W_b(z)$ est sans dimension, donc $T_{A\to B}$ a des unités de temps (comme requis).

Guide pratique

Entrées dont vous avez besoin.
- $T_{B\to C}$ : la demi-vie terminale du métabolite actif $B$ (à partir des données concentration–temps).
- $t_\text{peak}$ : le temps de l’effet maximal (ou du $B$ maximal si disponible). Si l’effet est utilisé, assurez-vous que l’effet est directement proportionnel à $B$ sans délai (pas d’hystérésis).
Calculez le $a$ sans dimension.
\[ a=\frac{(\ln 2)\,t_\text{peak}}{T_{B\to C}}. \]
Choisissez la branche.
- Si $a>1$ : utilisez $W_0$ (formation plus lente que l’élimination).
- Si $a<1$ : utilisez $W_{-1}$ (formation plus rapide que l’élimination).
- Si $a=1$ : $T_{A\to B}=T_{B\to C}$.
Évaluez $T_{A\to B}$.
- Calculez $z=-a e^{-a}$, puis évaluez $W_b(z)$ dans votre bibliothèque mathématique ; enfin calculez $T_{A\to B}=-(\ln 2)\,t_\text{peak}/W_b(z)$.
Vérifications de cohérence.
- Si $t_\text{peak}$ est très tôt par rapport à $T_{B\to C}$ (petit $a$), vous devriez obtenir $T_{A\to B}\ll T_{B\to C}$.
- Si $t_\text{peak}$ est très tard, vous devriez obtenir $T_{A\to B}\gg T_{B\to C}$.

Cas limites et hypothèses

Taux égaux $k_f=k_e$ : traité comme ci-dessus avec $t_\text{peak}=1/k_e$.
Portée du modèle. La dérivation suppose :
- Une cinétique du premier ordre (linéaire) pour la formation et l’élimination.
- Un seul compartiment bien mélangé pour $B$ (pas de délais distributionnels).
- L’effet PD est proportionnel à $B$ sans délai (pas de compartiment d’effet / tolérance). Les violations (par exemple, métabolisme saturable, distribution multi-compartiments, réponse indirecte) décaleront $t_\text{peak}$ et briseront l’application.
Bruit de mesure. $t_\text{peak}$ estimé à partir d’un échantillonnage clairsemé peut être biaisé. Lissez ou ajustez par modèle l’évolution temporelle de $B$ (ou de l’effet) pour estimer le vrai temps de pic.

Code Python pour calculer la demi-vie d’absorption

compute_formation_half_life.py

#!/usr/bin/env python3
# Calcule T_{A->B} à partir de T_{B->C} et t_peak en utilisant Lambert W (NumPy + SciPy).
import numpy as np
from scipy.special import lambertw

def formation_half_life(T_BC, t_peak):
    """
    Retourne T_{A->B} (demi-vie de formation) dans la chaîne A->B->C, étant donné :
      - T_BC   : demi-vie de B->C (mêmes unités de temps que la sortie souhaitée)
      - t_peak : temps du pic de B (et du pic d'effet si E ∝ B)

    Forme analytique :
        T_AB = - (ln 2) * t_peak / W_b( -a * e^{-a} ),
      où a = (ln 2) * t_peak / T_BC
    Règle de branche (solution réelle) :
        si a > 1 -> utiliser W_0   (formation plus lente que l'élimination)
        si a < 1 -> utiliser W_{-1} (formation plus rapide que l'élimination)
        si a == 1 -> T_AB = T_BC
    """
    T_BC_arr, t_peak_arr = np.broadcast_arrays(np.asarray(T_BC, float),
                                               np.asarray(t_peak, float))

    if np.any(T_BC_arr <= 0) or np.any(t_peak_arr <= 0):
        raise ValueError("T_BC et t_peak doivent être positifs.")

    ln2 = np.log(2.0)
    a = ln2 * t_peak_arr / T_BC_arr
    z = -a * np.exp(-a)

    T_AB = np.empty_like(a, dtype=float)

    # Cas spécial : a == 1 -> T_AB == T_BC
    mask_eq = np.isclose(a, 1.0, rtol=1e-12, atol=0.0)
    if np.any(mask_eq):
        T_AB[mask_eq] = T_BC_arr[mask_eq]

    # a > 1 -> branche principale W_0 -> k_f < k_e -> T_AB > T_BC
    mask_gt = a > 1.0
    if np.any(mask_gt):
        W0 = lambertw(z[mask_gt], 0)         # type complexe
        T_AB[mask_gt] = -(ln2 * t_peak_arr[mask_gt]) / W0.real

    # a < 1 -> branche inférieure W_{-1} -> k_f > k_e -> T_AB < T_BC
    mask_lt = a < 1.0
    if np.any(mask_lt):
        Wm1 = lambertw(z[mask_lt], -1)       # type complexe
        T_AB[mask_lt] = -(ln2 * t_peak_arr[mask_lt]) / Wm1.real

    # Retourne un scalaire si les entrées étaient des scalaires
    if np.isscalar(T_BC) and np.isscalar(t_peak):
        return float(T_AB.item())
    return T_AB


def formation_rate_constant(T_BC, t_peak):
    """Commodité : retourne k_f = ln(2) / T_{A->B} (même comportement de diffusion)."""
    T_AB = formation_half_life(T_BC, t_peak)
    return np.log(2.0) / T_AB


if __name__ == "__main__":
    # Exemples
    # 1) a == 1 -> T_AB == T_BC
    T_BC = 10.0
    t_peak = T_BC / np.log(2.0)
    print("Exemple a=1 :", formation_half_life(T_BC, t_peak))

    # 2) formation plus lente (a > 1) -> T_AB > T_BC
    print("Exemple a>1 :", formation_half_life(8.0, 20.0))

    # 3) formation plus rapide (a < 1) -> T_AB < T_BC
    print("Exemple a<1 :", formation_half_life(12.0, 2.0))

    # 4) entrées vectorisées
    T_BC_vec = np.array([10.0, 8.0, 12.0])
    t_peak_vec = np.array([10.0/np.log(2.0), 20.0, 2.0])
    print("Vectorisé :", formation_half_life(T_BC_vec, t_peak_vec))

#!/usr/bin/env python3
# Calcule T_{A->B} à partir de T_{B->C} et t_peak en utilisant Lambert W (NumPy + SciPy).
import numpy as np
from scipy.special import lambertw

def formation_half_life(T_BC, t_peak):
    """
    Retourne T_{A->B} (demi-vie de formation) dans la chaîne A->B->C, étant donné :
      - T_BC   : demi-vie de B->C (mêmes unités de temps que la sortie souhaitée)
      - t_peak : temps du pic de B (et du pic d'effet si E ∝ B)

    Forme analytique :
        T_AB = - (ln 2) * t_peak / W_b( -a * e^{-a} ),
      où a = (ln 2) * t_peak / T_BC
    Règle de branche (solution réelle) :
        si a > 1 -> utiliser W_0   (formation plus lente que l'élimination)
        si a < 1 -> utiliser W_{-1} (formation plus rapide que l'élimination)
        si a == 1 -> T_AB = T_BC
    """
    T_BC_arr, t_peak_arr = np.broadcast_arrays(np.asarray(T_BC, float),
                                               np.asarray(t_peak, float))

    if np.any(T_BC_arr <= 0) or np.any(t_peak_arr <= 0):
        raise ValueError("T_BC et t_peak doivent être positifs.")

    ln2 = np.log(2.0)
    a = ln2 * t_peak_arr / T_BC_arr
    z = -a * np.exp(-a)

    T_AB = np.empty_like(a, dtype=float)

    # Cas spécial : a == 1 -> T_AB == T_BC
    mask_eq = np.isclose(a, 1.0, rtol=1e-12, atol=0.0)
    if np.any(mask_eq):
        T_AB[mask_eq] = T_BC_arr[mask_eq]

    # a > 1 -> branche principale W_0 -> k_f < k_e -> T_AB > T_BC
    mask_gt = a > 1.0
    if np.any(mask_gt):
        W0 = lambertw(z[mask_gt], 0)         # type complexe
        T_AB[mask_gt] = -(ln2 * t_peak_arr[mask_gt]) / W0.real

    # a < 1 -> branche inférieure W_{-1} -> k_f > k_e -> T_AB < T_BC
    mask_lt = a < 1.0
    if np.any(mask_lt):
        Wm1 = lambertw(z[mask_lt], -1)       # type complexe
        T_AB[mask_lt] = -(ln2 * t_peak_arr[mask_lt]) / Wm1.real

    # Retourne un scalaire si les entrées étaient des scalaires
    if np.isscalar(T_BC) and np.isscalar(t_peak):
        return float(T_AB.item())
    return T_AB


def formation_rate_constant(T_BC, t_peak):
    """Commodité : retourne k_f = ln(2) / T_{A->B} (même comportement de diffusion)."""
    T_AB = formation_half_life(T_BC, t_peak)
    return np.log(2.0) / T_AB


if __name__ == "__main__":
    # Exemples
    # 1) a == 1 -> T_AB == T_BC
    T_BC = 10.0
    t_peak = T_BC / np.log(2.0)
    print("Exemple a=1 :", formation_half_life(T_BC, t_peak))

    # 2) formation plus lente (a > 1) -> T_AB > T_BC
    print("Exemple a>1 :", formation_half_life(8.0, 20.0))

    # 3) formation plus rapide (a < 1) -> T_AB < T_BC
    print("Exemple a<1 :", formation_half_life(12.0, 2.0))

    # 4) entrées vectorisées
    T_BC_vec = np.array([10.0, 8.0, 12.0])
    t_peak_vec = np.array([10.0/np.log(2.0), 20.0, 2.0])
    print("Vectorisé :", formation_half_life(T_BC_vec, t_peak_vec))

Exemple d’utilisation : Absorption de la quétiapine

Pour cet exemple, nous utiliserons l’absorption de la quétiapine orale (métabolite A) en quétiapine (métabolite B).

Les données pharmacologiques expérimentales selon NCBI sont :

environ 1 à 2 heures => nous utiliserons la moyenne 1,5h pour $T_{A\to B}$
La demi-vie rapportée (t1/2) pour la quétiapine est d'environ 7 heures. pour $T_{B\to C}$

Comme le montre le graphique suivant, la concentration plasmatique maximale de la littérature correspond à la simulation analytique basée sur la demi-vie d’absorption calculée.

Quetiapine pharmacokinetics.svg

Sortie :

quetiapine_example_output.txt

T_BC = 6.0 h, Tmax = 1.5 h
Computed uptake half-life T_AB ≈ 0.3424 h (20.5 min)
kf = 2.0246 1/h, ke = 0.1155 1/h

T_BC = 6.0 h, Tmax = 1.5 h
Computed uptake half-life T_AB ≈ 0.3424 h (20.5 min)
kf = 2.0246 1/h, ke = 0.1155 1/h

plot_quetiapine_pk.py

#!/usr/bin/env python3
import numpy as np, math
import matplotlib.pyplot as plt
plt.style.use("ggplot")

from scipy.special import lambertw

ln2 = math.log(2.0)
T_BC_hours = 6.0    # demi-vie d'élimination de la quétiapine (B -> élimination)
t_peak_hours = 1.5  # Tmax (temps jusqu'à la concentration plasmatique maximale)

def compute_T_AB(T_BC, t_peak):
    """Retourne la demi-vie d'absorption T_AB à partir de T_BC et t_peak.
       Utilise Lambert W si SciPy est disponible ; sinon bissection sur a=ln r/(r-1)."""
    if T_BC <= 0 or t_peak <= 0:
        raise ValueError("T_BC et t_peak doivent être positifs.")
    a = ln2 * t_peak / T_BC
    if abs(a - 1.0) <= 1e-12:
        return T_BC  # cas des taux égaux

    # Calcule en utilisant Lambert W (SciPy présent)
    z = -a * math.exp(-a)
    branch = 0 if a > 1 else -1
    W = lambertw(z, branch)
    return float(-(ln2 * t_peak) / W.real)

# Calcule les demi-vies et les taux
T_AB_hours = compute_T_AB(T_BC_hours, t_peak_hours)
ke = ln2 / T_BC_hours
kf = ln2 / T_AB_hours

print(f"T_BC = {T_BC_hours} h, Tmax = {t_peak_hours} h")
print(f"Demi-vie d'absorption calculée T_AB ≈ {T_AB_hours:.4f} h ({T_AB_hours*60:.1f} min)")
print(f"kf = {kf:.4f} 1/h, ke = {ke:.4f} 1/h")

# Grille temporelle et quantités normalisées (A(0)=1)
t = np.linspace(0.0, 24.0, 1000)  # heures
A = np.exp(-kf * t)
B = (kf / (ke - kf)) * (np.exp(-kf * t) - np.exp(-ke * t))  # forme Bateman normalisée

# Tracé
plt.figure(figsize=(8,5))
plt.plot(t, A, label="A(t) : quétiapine orale (normalisée)")
plt.plot(t, B, label="B(t) : quétiapine (normalisée)")
plt.axvline(t_peak_hours, linestyle="--", label=f"Tmax de la littérature ≈ {t_peak_hours:.2f} h")
plt.xlabel("Temps (heures)")
plt.ylabel("Quantité / concentration normalisée (u.a.)")
plt.title("Quétiapine — demi-vie d'absorption dérivée et courbes PK normalisées")
plt.legend()
plt.tight_layout()
plt.savefig("Quetiapine pharmacokinetics.svg")

#!/usr/bin/env python3
import numpy as np, math
import matplotlib.pyplot as plt
plt.style.use("ggplot")

from scipy.special import lambertw

ln2 = math.log(2.0)
T_BC_hours = 6.0    # demi-vie d'élimination de la quétiapine (B -> élimination)
t_peak_hours = 1.5  # Tmax (temps jusqu'à la concentration plasmatique maximale)

def compute_T_AB(T_BC, t_peak):
    """Retourne la demi-vie d'absorption T_AB à partir de T_BC et t_peak.
       Utilise Lambert W si SciPy est disponible ; sinon bissection sur a=ln r/(r-1)."""
    if T_BC <= 0 or t_peak <= 0:
        raise ValueError("T_BC et t_peak doivent être positifs.")
    a = ln2 * t_peak / T_BC
    if abs(a - 1.0) <= 1e-12:
        return T_BC  # cas des taux égaux

    # Calcule en utilisant Lambert W (SciPy présent)
    z = -a * math.exp(-a)
    branch = 0 if a > 1 else -1
    W = lambertw(z, branch)
    return float(-(ln2 * t_peak) / W.real)

# Calcule les demi-vies et les taux
T_AB_hours = compute_T_AB(T_BC_hours, t_peak_hours)
ke = ln2 / T_BC_hours
kf = ln2 / T_AB_hours

print(f"T_BC = {T_BC_hours} h, Tmax = {t_peak_hours} h")
print(f"Demi-vie d'absorption calculée T_AB ≈ {T_AB_hours:.4f} h ({T_AB_hours*60:.1f} min)")
print(f"kf = {kf:.4f} 1/h, ke = {ke:.4f} 1/h")

# Grille temporelle et quantités normalisées (A(0)=1)
t = np.linspace(0.0, 24.0, 1000)  # heures
A = np.exp(-kf * t)
B = (kf / (ke - kf)) * (np.exp(-kf * t) - np.exp(-ke * t))  # forme Bateman normalisée

# Tracé
plt.figure(figsize=(8,5))
plt.plot(t, A, label="A(t) : quétiapine orale (normalisée)")
plt.plot(t, B, label="B(t) : quétiapine (normalisée)")
plt.axvline(t_peak_hours, linestyle="--", label=f"Tmax de la littérature ≈ {t_peak_hours:.2f} h")
plt.xlabel("Temps (heures)")
plt.ylabel("Quantité / concentration normalisée (u.a.)")
plt.title("Quétiapine — demi-vie d'absorption dérivée et courbes PK normalisées")
plt.legend()
plt.tight_layout()
plt.savefig("Quetiapine pharmacokinetics.svg")

Check out similar posts by category: Bioinformatics, Pharmacokinetics, Mathematics

If this post helped you, please consider buying me a coffee or donating via PayPal to support research & publishing of new posts on TechOverflow

Buy me a coffee

Pharmacocinétique : Une solution analytique pour calculer la demi-vie d'absorption à partir du temps de concentration plasmatique maximale