Învățarea simularii PK/PD: o analiză Monte Carlo pentru începători cu mrgsolve în R

🧪 Scufundarea în PK/PD pentru prima dată – simularea ceftriaxonei cu mrgsolve în R. Nivelurile de medicamente gratuite au fost… surprinzător de ridicate? Chiar și l-am împins la dozarea q48h din curiozitate, iar rezultatele m-au lăsat cu mai multe întrebări decât răspunsuri 🤔📈

Motivații

Învățarea farmacocineticii (PK) și a farmacodinamicii (PD) a fost întotdeauna un interes pentru mine. Este întotdeauna o provocare să citim aceste documente PK populaționale cu toate numerele etc. Care este o modalitate mai bună de a ne scufunda în suprafața acestora și de a vedea dacă măcar putem ști cum să codificăm o simulare pentru a obține probabilitatea atingerii țintei (PTA) a diferitelor concentrații minime inhibitorii (mic) și să învățăm elementele de bază prin cod! Să mergem mai departe!

Disclaimer:

Nu sunt farmacist și nici expert în PK/PD. Aceasta este o documentație pentru propria mea învățare și numai în scopuri educaționale. Nu un sfat medical. Dacă ați observat ceva în neregulă aici, vă rog să-mi spuneți!

Obiective:

Ce este populația PK

Farmacocinetica populației (popPK) este o abordare statistică care descrie modul în care medicamentele se comportă în organism pe grupuri de oameni, ținând cont de variabilitatea dintre indivizi. În loc să studieze o singură persoană în mod intensiv, popPK analizează date rare de la mulți pacienți pentru a înțelege comportamentul tipic la medicamente și de ce oamenii diferă în ceea ce privește expunerea la medicamente.

Vom folosi Garot D și colab. Farmacocinetica populației a ceftriaxonei la pacienții septici critici: o reevaluare ca exemplu pentru învățare.

Care sunt parametrii de interes pe o lucrare?

Din hârtie, putem vedea că există o mulțime de parametri și numere. Dar care sunt parametrii de interes? Ne vom concentra pe următorii parametri:

Privind la lor table 3vedem aceste valori:
CL = (theta_1 + theta_2 . (CL_{cr}/4.26))
$theta_1$ : Componentă de clearance non-renală (de bază).
$theta_2$ : Coeficientul de scalare a clearance-ului renal
V1 : Volumul de distribuție al compartimentului central.
V2 : Volumul de distribuţie al compartimentului periferic.
Q : Degajare intercompartimentală.
$omega^2 (CL)$ : Variabilitatea clearance-ului între subiecte.
$omega^2 (V1)$ : Variabilitatea între subiecte a volumului de distribuție a compartimentului central.
$omega^2 (V2)$ : Variabilitatea între subiecți a volumului de distribuție a compartimentului periferic.

Aceștia sunt parametrii pe care îi vom folosi în modelul nostru mrgsolve. Întotdeauna m-am întrebat ce reprezintă acești parametri și a fost puțin dificil de conceptualizat până când ne-am aruncat în cod și în sfârșit am înțeles rațiunea! Este un model cu efect mixt în care estimările au fost modelate în funcție de efectul fix (theta) și efectul aleator (eta), veți vedea acest lucru în cod mai târziu. Efectul fix reprezintă valoarea tipică a parametrului în populație, în timp ce efectul aleator reprezintă variabilitatea dintre indivizi. Se presupune că efectul aleatoriu este distribuit în mod normal cu o medie de zero și o varianță de omega pătrat.

Dacă ar fi să desenăm o diagramă a celor de mai sus, va arăta cam așa:

Unul începe cu administrarea medicamentelor în compartimentul central, apoi merge fie în compartimentul periferic (țesut etc.) și în clearance-ul medicamentelor. Observați că Q este un flux bidirecțional între central și periferic, în timp ce toate celelalte direcții sunt fie în centru, fie în afară de la centru până la clearance. Acest lucru este foarte util pentru mine pentru a obține o înțelegere superficială a distribuției. Să continuăm cu codul!

Să codificăm

library(mrgsolve)
library(tidyverse)

mod <- mcode(model = "ceftriaxone", code='
$PARAM
theta1 = 0.56,
theta2 = 0.32,
CLcr   = 4.26,   // median creatinine clearance of 68.5 ml min-1, hence ~4.26 L hr-1
V1     = 10.3,
V2     = 7.35,
Q      = 5.28,  
fu     = 0.10  // fraction of unbound, picked a static value from package insert range

$CMT CENT PERI

$MAIN
double CL  = theta1 + theta2 * (CLcr / 4.26);  // technically we could use 0.88 as reported on their result section
double CLi = CL * exp(ETA(1));                  // ETA here means log normal distibution of mean 
double V1i = V1 * exp(ETA(2));     
double V2i = V2 * exp(ETA(3));

$OMEGA
0.24   // omega2(CL) from table
0.23
0.42

$SIGMA
0.0576   // √0.0576 = 0.24, 

$ODE
dxdt_CENT = -(CLi/V1i)*CENT - (Q/V1i)*CENT + (Q/V2i)*PERI;
dxdt_PERI =  (Q/V1i)*CENT   - (Q/V2i)*PERI;

$TABLE
double Cp_total = (CENT / V1i)*(1+EPS(1)); 
double Cp_free = fu * Cp_total;

$CAPTURE Cp_free
')

dosing <- ev(amt = 2000, rate = 4000, ii = 24, addl = 2, cmt = "CENT")

set.seed(1)
sims <- mod |>
  ev(dosing) |>
  mrgsim(nid = 1000, end = 72, delta = 0.25) |>
  as_tibble()

ETA în cele de mai sus înseamnă efectul aleatoriu, care se presupune a fi distribuit în mod normal cu o medie de zero și o varianță de omega pătrat. ETA este o literă greacă (eh-ta). EPS aici este Epsilon.

Pentru evcantitatea este dozată în mg; rate este suma dată pe oră; ii este frecvența; addl este numărul de doze suplimentare; cmt este compartimentul în care se administrează doza. În acest caz, administrăm 2000 mg de ceftriaxonă sub formă de perfuzie de 30 de minute la fiecare 24 de ore pentru 3 doze (1 doză inițială + 2 doze suplimentare) în compartimentul central.

Apoi, va trebui să setăm sămânța pentru reproductibilitate, apoi în modelul dumneavoastră inițial cu dozare nid este numărul de persoane pe care doriți să le simulați, end este ora de încheiere a simulării în ore și delta este intervalul de timp pentru ieșirea de simulare în ore. În acest caz, simulăm 1000 de indivizi timp de 72 de ore cu un interval de timp de 0,25 ore (15 minute).

În continuare, vom calcula probabilitatea de atingere a țintei (PTA) pentru diferite valori ale concentrației minime inhibitorii (MIC). PTA este probabilitatea ca concentrația de medicament liber să depășească CMI pentru un anumit procent din intervalul de dozare.

MIC <- 1

print(paste0("Probability of Target Attainment: ", sims |>
  filter(time >= 48) |>
  group_by(ID) |>
  summarise(fT = mean(Cp_free > MIC)) |>
  summarise(PTA = mean(fT >= 0.50)) |>
  pull()))

## (1) "Probability of Target Attainment: 0.996"

sims |>
  ggplot(aes(x=time,y=Cp_free,group=ID)) +
  geom_line(alpha=0.01) +
  geom_hline(yintercept = MIC, color = "red") +
  theme_bw()

În cele de mai sus, alegem microfon de 1, timp de oprire filtrat după 48 de ore pentru starea de echilibru, apoi calculăm ceftriaxona liberă medie care este deasupra microfonului, apoi evaluăm media timpilor în care ceftriaxona liberă este peste 50% per subiect simulat. Putem observa că probabilitatea de atingere a țintei este de aproximativ 99,6%. De asemenea, putem vizualiza concentrația de ceftriaxonă liberă în timp cu o linie roșie care indică microfonul de 1. Nu prea ponosit! Acum să evaluăm când există o diferență între CrCl și albumină. Vă scutesc de cod.

Modificări ale CrCl

library(glue)

crcl_vec <- c(1.8,7.2,10.8)
crcl_vec_i <- c(30, 120, 180)

# 30 ml/min = 1.8 L/hr
# 120 ml/min = 7.2 L/hr
# 180 ml/min = 10.8 L/hr

for (crcl in crcl_vec) {

mod <- mcode(model = "ceftriaxone", code=glue('
$PARAM
theta1 = 0.56,
theta2 = 0.32,
CLcr   = {crcl},
V1     = 10.3,
V2     = 7.35,
Q      = 5.28,  
fu     = 0.10  // fraction of unbound, picked a static value from package insert range

$CMT CENT PERI

$MAIN
double CL  = theta1 + theta2 * (CLcr / 4.26);  // technically we could use 0.88 as reported on their result section
double CLi = CL * exp(ETA(1));                  // ETA here means log normal distibution of mean 
double V1i = V1 * exp(ETA(2));     
double V2i = V2 * exp(ETA(3));

$OMEGA
0.24   // omega2(CL) from table
0.23
0.42

$SIGMA
0.0576   // √0.0576 = 0.24, 

$ODE
dxdt_CENT = -(CLi/V1i)*CENT - (Q/V1i)*CENT + (Q/V2i)*PERI;
dxdt_PERI =  (Q/V1i)*CENT   - (Q/V2i)*PERI;

$TABLE
double Cp_total = (CENT / V1i)*(1+EPS(1)); 
double Cp_free = fu * Cp_total;

$CAPTURE Cp_free
',crcl))

dosing <- ev(amt = 2000, rate = 4000, ii = 24, addl = 2, cmt = "CENT")

set.seed(1)
sims <- mod |>
  ev(dosing) |>
  mrgsim(nid = 1000, end = 72, delta = 0.25) |>
  as_tibble()

MIC <- 1

pta <- paste0("Probability of Target Attainment: ", sims |>
  filter(time >= 48) |>
  group_by(ID) |>
  summarise(fT = mean(Cp_free > MIC)) |>
  summarise(PTA = mean(fT >= 0.50)) |>
  pull(), " ,CrCl: ", crcl_vec_i(crcl_vec==crcl), "ml/min")


plot <- sims |>
  ggplot(aes(x=time,y=Cp_free,group=ID)) +
  geom_line(alpha=0.01) +
  geom_hline(yintercept = MIC, color = "red") +
  theme_bw() +
  ggtitle(pta)

plot(plot)
}

E interesant! Acest lucru are sens, creșterea CrCl va crește clearance-ul ceftriaxonei, prin urmare scăderea PTA. Totusi e destul de bine! Totuși, ce este considerat acceptabil? 90%? 70%? 50%? De asemenea, PTA de mai sus se bazează pe 50% din timp ceftriaxona liberă este peste MIC. Atunci care este numărul acceptabil pentru asta? 🤷‍♂️ Ce se întâmplă dacă, din moment ce ceftriaxona este legată de albumină, dacă modelăm și albumină în model?

Modificări cu hipoalbuminemia

Observați că modelul nostru inițial a avut o fracțiune fixă ​​de nelegat (fu) de 0,1, care este mijlocul intervalului raportat în prospect. Cu toate acestea, la pacienții în stare critică, hipoalbuminemia este frecventă și poate duce la o creștere a fracției de medicament nelegat, ceea ce poate afecta farmacocinetica și farmacodinamia ceftriaxonei. Să vedem cum putem modela acest lucru în codul nostru mrgsolve. Din lucrarea din secțiunea metode, au folosit formula de mai jos pentru a estima ceftriaxona liberă din ceftriaxona totală:

Vom adăuga asta la modelul nostru și vom ajusta np (concentrația totală de situsuri de legare a proteinelor) conform estimării cu albumină mai mică (np = 295), acest număr a fost din nou din lucrarea din porțiunea de discuție în care albumina lor mediană a fost de ~ 25 g/L.

mod <- mcode(model = "ceftriaxone", code='
$PARAM
theta1 = 0.56,
theta2 = 0.32,
CLcr   = 4.26,
V1     = 10.3,
V2     = 7.35,
Q      = 5.28,
np     = 517,  
kaff   = 0.0367 

$OMEGA
0.24
0.23
0.42

$SIGMA
0.0576

$CMT CENT PERI

$GLOBAL
double solveFree(double CTOT, double np, double kaff) {
  double cf   = (-(np+1/kaff-CTOT)+sqrt(pow(np+1/kaff-CTOT,2.0)+(4.0*CTOT/kaff)));
  return cf > 0 ? cf : 0;
}

$MAIN
double CL  = theta1 + theta2 * (CLcr / 4.26);
double CLi = CL * exp(ETA(1));
double V1i = V1 * exp(ETA(2));
double V2i = V2 * exp(ETA(3));

$ODE
double CTOT  = CENT / V1i;           // renamed: avoid clash with $TABLE
double CFREE = solveFree(CTOT, np, kaff);

dxdt_CENT = -CLi * CFREE
            - (Q / V1i) * CENT
            + (Q / V2i) * PERI;

dxdt_PERI =  (Q / V1i) * CENT
            - (Q / V2i) * PERI;

$TABLE
double CTOTAL      = CENT / V1i;      // notice this is not CTOT
double Cp_free     = solveFree(CTOTAL, np, kaff);
double Cp_bound    = CTOTAL - Cp_free;
double FU          = Cp_free / (CTOTAL + 1e-9);
double Cp_obs      = CTOTAL * (1 + EPS(1));

$CAPTURE CTOTAL Cp_free Cp_bound FU Cp_obs
')

dosing <- ev(amt = 2000, rate = 4000, ii = 24, addl = 2, cmt = "CENT")

set.seed(1)
sims <- mod |>
  param(np = 295) |>
  ev(dosing) |>
  mrgsim(nid = 1000, end = 72, delta = 0.25) |>
  as_tibble()

MIC <- 1

pta <- paste0("Probability of Target Attainment: ", sims |>
  filter(time >= 48) |>
  group_by(ID) |>
  summarise(fT = mean(Cp_free > MIC)) |>
  summarise(PTA = mean(fT >= 0.50)) |>
  pull(),", Albumin: ~25g/L, CrCl: ~63 ml/min")

plot <- sims |>
  ggplot(aes(x=time,y=Cp_free,group=ID)) +
  geom_line(alpha=0.01) +
  geom_hline(yintercept = MIC, color = "red") +
  # geom_text(aes(x=20,y=150,label=pta)) +
  theme_bw() +
  ggtitle(pta)

plot(plot)

Wow, e interesant! După ce ne-am potrivit corect în estimarea ceftriaxonei libere, aceasta a îmbunătățit de fapt PTA chiar și atunci când albumina este mai mică. Ce se întâmplă dacă facem albumina chiar mai mică la ~15g/L (np=~172) și creștem CrCl la 180 ml/min și creștem fT >= 0,7 (mai mult de 70% din timp ceftriaxona liberă este peste mic) și vedem dacă vom putea curăța medicamentul mai repede?

PTA este încă 100% !?!?! wow, ceftriaxona 2g chiar este o bestie! Hmmm.. Ceftriaxona liberă este FOARTE mare, în jur de ~200-300, putem simula o dozare la 48 de ore și să vedem cum este PTA, chiar și pentru scenariul nostru cel mai rău, albumină scăzută, CrCl ridicat și acoperire ft>mic >= 70%?

?q48 Dozare

dosing <- ev(amt = 2000, rate = 4000, ii = 48, addl = 2, cmt = "CENT")

set.seed(1)
sims <- mod |>
  param(np = 172, CLcr = 10.8) |>
  ev(dosing) |>
  mrgsim(nid = 1000, end = 144, delta = 0.25) |>
  as_tibble()

MIC <- 1

pta <- paste0("PTA: ", sims |>
  filter(time >= 48) |>
  group_by(ID) |>
  summarise(fT = mean(Cp_free > MIC)) |>
  summarise(PTA = mean(fT >= 0.70)) |>
  pull(),", Albumin: ~15g/L, CrCl: ~180 ml/min, fT > mic >= 70%, q48h dosing")

plot <- sims |>
  ggplot(aes(x=time,y=Cp_free,group=ID)) +
  geom_line(alpha=0.01) +
  geom_hline(yintercept = MIC, color = "red") +
  theme_bw() +
  ggtitle(pta)

plot(plot)

Serios!? PTA este încă atât de mare!? Ce înseamnă asta de fapt? Există literatură despre asta? Poate că codul meu nu este corect… 🤔🤷‍♂️

Să ne uităm la fT > mic >= 99%.

## (1) "PTA: 0.979, Albumin: ~15g/L, CrCl: ~180 ml/min, fT > mic >= 99%, q48h dosing"

Dacă știi ceva despre asta, te rog să-mi spui! aceasta este pentru organismul cu mic <= 1, ceftriaxonă 2g. Din nou, rețineți că acest lucru este doar în scopuri educaționale și de învățare. Descoperirea pe care am primit-o mai sus este doar o explorare curioasă. Mă întreb dacă există o eroare de codare din partea mea. Faceți clic pe code mai sus pentru a extinde pentru detalii. De asemenea, mă întreb dacă majoritatea încercărilor pe care le-am avut înainte s-au bazat pe un microfon mai mare, în timp ce microfonul din zilele noastre pentru ceftriaxonă este în principal <= 1. 🤔

Oportunități de îmbunătățire

• Aflați cum modelează popPK, acest lucru ne va ajuta cu adevărat să înțelegem cum au obținut acele estimări theta și omega
• Nu prea înțeleg încă partea sigma, mă voi scufunda în asta data viitoare, mai ales când estim aceste valori
• Încercați să aflați alte proprietăți, cum ar fi AUC/mic, Cmax/mic etc. și vedeți cum se modifică PTA
• Aflați mai multe din literatura de specialitate care este preferată în ceea ce privește nivelul acceptabil de medicație gratuită peste microfon și PTA acceptabil

Lecții învățate

• am învățat codificarea modelului mrgsolve (folosește cpp)
• am învățat câteva ecuații de bază pk/pd, popPK
• a aflat despre cele 2 compartimente
• am găsit un rezultat neașteptat pentru dozarea q48 prin simulare, încă nu sunt sigur dacă acesta este ceva real/adevărat
• am învățat că theta nu este legat de central/periferic, mai degrabă theta1 este clearance-ul de bază și theta2 este? scalarea renală

Daca iti place acest articol: