Jurnal Farmasi Galenika ( Galenika Journal of Pharmacy ) 2018; 4 (1): 26 – 35 ISSN : 2442-8744 (electronic) http://jurnal.untad.ac.id/jurnal/index.php/Galenika/index DOI : 10.22487/j24428744. 2018.v4.i1.9777

## Pemodelan  Farmakokinetika  Berbasis  Populasi  dengan  R:  Model  Dua Kompartemen Ekstravaskuler
( Population-Based  Pharmacokinetics  Modeling  with  R:  Two  Compartment  Extravascular Model )
Dion Notario *
Program Studi Farmasi, Fakultas Kedokteran, Universitas Katolik Indonesia Atma Jaya, Jakarta, Indonesia
Article Info: Received: 16 Maret 2018 in revised form: 30 Maret 2018 Accepted: 2 April 2018
Available Online: 2 April 2018

## ABSTRACT
A Tutorial of two-compartment extravascular population-based pharmacokinetics  modeling  was  performed  by  differential  equations  and non-linear  mixed  effect  model  approach.  First,  three-level  differential equations  of  two-compartment  pharmacokinetics  were  generated.  Then, covariate and non-covariate models were developed by nlmeODE and nlme packages installed in R. The best model was selected according to AIC, BIC, and  LogLik  value.  A  model  without  covariates  model  was  selected  as  the best model. The selected model showed a goodness of fit with experimental dataset  and  residual  plot of  the  model  revealed  that  no  violations  of  model assumtions.    In  conclusion,  nlme  and  nlmeODE  is  capable  to  generate  an adequate predictive model of two-compartment population-based pharmacokinetics for extravascular route.
Keywords :
Population-based pharmacokinetics Two-compartement Differential equation nlme and nlmeODE

## Corresponding Author:
Dion Notario Program Studi Farmasi, FK, Universitas Katolik Indonesia Atma Jaya Jakarta Utara,14440, Jakarta, Indonesia Mobile : 081804175860 Email: dion.natitivate@gmail.com
Copyright © 2017 JFG-UNTAD
This open access article is distributed under a Creative Commons Attribution (CC-BY-NC-SA) 4.0 International license.

## How to cite (APA 6th Style):
Notario,  D.  (2018).  Pemodelan  Farmakokinetika  Berbasis  Populasi  dengan  R:  Model  Dua  Kompartemen Ekstravaskuler. Jurnal Farmasi Galenika : Galenika Journal of Pharmacy, 4(1), 26-35.
doi:10.22487/j24428744.2017.v4.i1.9777

## ABSTRAK
Sebuah  tutorial  farmakokinetika  dua  kompartemen  ekstravaskuler  berbasis  populasi  dilakukan dengan  pendekatan  persamaan  diferensial  dan  model  efek  campuran  non  linear.  Pertama,  disusun persamaan  diferensial  tiga  tingkat  dari  model  dua  kompartemen  ekstravaskuler.Kemudian, dikembangkan  model  kovariat  dan  non-kovariat  dengan  paket  nlmeODE  dan  nlme  yang  terpasang dalam  perangkat  lunak  R.  Model  terbaik  dipilih  berdasarkan  nilai  AIC,  BIC,  dan  LogLik.  Sebuah model  tanpa  kovariat  dipilih  sebagai  model  terbaik.  Model  yang  terpilih  menunjukkan  kesesuaian fungsi  yang  baik  dengan  data  eksperimental  dan  plot  residual  model  menunjukkan  bahwa  tidak terdapat  satu  pun  asumsi  yang  dilanggar.  Dengan  demikian,  dapat  disimpulkan  bahwa  nlme  dan nlmeODE mampu menghasilkan model prediktif farmakokinetika dua kompartemen berbasis populasi yang memadai untuk rute ekstravaskuler.
Kata  Kunci  :  Farmakokinetika  berbasis  populasi,  Dua  Kompartemen,    Persamaan  differensial,  nlme dan nlmeODE.

## PENDAHULUAN
Farmakokinetika adalah ilmu yang mempelajari  secara  khusus  perubahan  jumlah  obat dalam  tubuh  sebagai  fungsi  waktu  (Jambhekar  & Breen,  2012;  Shargel et  al ,  2012).  Dengan  kata  lain, dalam  pokok  bahasan  farmakokinetika  dilakukan kajian-kajian  terhadap  fenomena  absorbsi,  distribusi, dan eliminasi obat secara kuantitatif. Oleh karena itu, dalam penelitian-penelitian farmakokinetika dikembangkan berbagai macam model-model matematika  untuk  menjelaskan  proses  perjalanan obat di dalam tubuh.
Salah  satu  pemodelan  matematika  yang paling  umum  digunakan  untuk  mengkaji  profil farmakokinetika  adalah  model  kompartementeral. Beberapa  model  farmakokinetika  kompartementeral antara  lain  model  kompartemen  tunggal  dan  multi kompartemen  telah  dikenal  secara  luas.  Diantara ketiga  model  kompartemen  tersebut,  model  dua kompartemen mengalami perkembangan yang paling pesat  yang  ditunjukkan  dengan  peningkatan  jumlah publikasi yang paling tinggi sejak tahun 1970 sampai 2017 (Kovalchik,  2017;  R  Core  Team,  2015) .  Oleh karena  itu,  pengetahuan  yang  lebih  dalam  mengenai model  dua  kompartemen  diperlukan  bagi  seorang farmasis  antara  lain  untuk  melakukan  pemodelan farmakokinetika  maupun  pengaturan  dosis  obat  di bidang farmasi klinis.
Beberapa model matematik untuk melakukan pemodelan  farmakokinetika  dua  kompartemen  telah dipublikasikan  antara  lain  oleh  Boxenbaum et  al , (1974)  dan  Wagner  (1975).  Beberapa  perangkat lunak juga telah dikembangkan untuk mempermudah
proses  kalkulasi  parameter  farmakokinetika  antara lain  dengan  WinSAAM,  R,  dan  PKSolver  (Lee  & Lee,  2017;  Stefanovski et  al ,  2003;  Tornøe et  al , 2004; Zhang et al , 2010).
Meskipun demikian, dari paket-paket perangkat  lunak  yang  telah  dikembangkan,  terutama perangkat  lunak  sumber  terbuka,  belum  pernah dilakukan pembahasan secara komprehensif mengenai proses kalkulasi parameter farmakokinetika populasi untuk model dua kompartemen  ekstravaskuler  dengan  persamaan diferensial. Padahal dengan memanfaatkan persamaan  diferensial,  dapat  dikembangkan  berbagai model  farmakokinetika  secara  lebih  leluasa,  tanpa melakukan  penyelesaian  dengan  teknik  kalkulus integral.  Di  lain  pihak,  tren  farmakokinetika  saat  ini telah mengalami pergeseran dari teknik
farmakokinetika
konvensional menuju ke farmakokinetika  berbasis  populasi,  bahkan  FDA ( Food  and  Drug  Administration )  mewajibkan penggunaan farmakokinetika berbasis populasi dalam uji klinik fase 1 dan akhir fase 2b (FDA, 1999). Oleh karena itu, di dalam tulisan ini, dilakukan kajian yang berpusat  pada  penyelesaian  persamaan  diferensial untuk  pemodelan  farmakokinetika  dua  kompartemen berbasis  populasi  dengan  perangkat  lunak  sumber terbuka  yaitu  R  dengan  paket  nlme  dan  nlmeODE (Pinheiro & Bates, 2000; Tornøe et al ., 2004).

## PEMODELAN DUA KOMPARTEMEN DENGAN PERSAMAAN DIFERENSIAL
Obat-obat yang diberikan melalui rute injeksi intravena bolus dan mengikuti model dua
kompartemen  tidak  mengalami  proses  absorpsi melainkan masuk secara langsung ke dalam sirkulasi sistemik.  Selanjutnya,  obat  akan  mengalami  proses eliminasi  dan  distribusi  secara  bersamaan.  Pada  fase distribusi,  terjadi  proses  perpindahan  obat  dari kompartemen  pertama  yaitu  sirkulasi  sistemik menuju  ke  kompartemen  kedua  atau  jaringan.  Pada tahap ini proses distribusi lebih dominan dibandingkan  dengan  eliminasi  yang  mengakibatkan terjadinya  pembelokan  kurva  profil  farmakokinetika pada  fase  distribusi  (Gambar  1  intravena  A).  Ketika obat  memasuki  fase  eliminasi  (Gambar  1  intravena B),  obat  telah  mencapai  kesetimbangan  dengan kompartemen  kedua  (jaringan)  sehingga  proses eliminasi  menjadi  lebih  dominan  dibandingkan distribusi,  bahkan  kecepatan  distribusi  sudah  dapat diabaikan.  Oleh  karena  itu,  pada  fase  eliminasi nampak bahwa titik-titik data membentuk suatu garis lurus menurun pada kertas semilog yang menandakan bahwa proses tersebut hanya dipengaruhi oleh proses eliminasi obat yang mengikuti kinetika orde pertama.
Kinetika  perjalanan  obat  dalam  tubuh  pada obat-obat  yang  diberikan  secara  intravena  bolus dengan  model  dua  kompartemen  dapat  disusun berdasarkan  ilustrasi  pada  Gambar  2A.  Apabila perpindahan  antar  kompartemen  dianggap  mengikuti kinetika  orde  pertama  dan  tanda  positif  atau  negatif secara  berurutan  menggambarkan  perjalanan  obat masuk dan keluar dari masing-masing  kompartemen, maka    dapat  dituliskan  seperangkat  persamaan diferensial  untuk  menggambarkan  farmakokinetika dua  kompartemen  pada  pemberian  obat  secara intravena  bolus  seperti  yang  dituliskan  pada persamaan (1) dan (2).
Gambar 1. Ilustrasi profil farmakokinetika dua kompartemen pada rute intravena dan ekstravaskuler pada kertas semilog. A = fase eliminasi, B = fase distribusi, C = fase absorpsi
k 21 k 10 k 12 1 2 sirkulasi sistemik jaringan V 1 A k 21 k 10 k 12 1 2 k a Tempat absorpsi sirkulasi sistemik jaringan V 1 TA B
Gambar 2. Ilustrasi model dua kompartemen intravena bolus (A) dan ekstravaskuler (B).
𝑑𝐶 1 𝑑𝑡 = 𝑘 21 𝐶 2 − 𝑘 12 𝐶 1 − 𝑘 10 𝐶 1 (1) 𝑑𝐶 2 𝑑𝑡 = 𝑘 12 𝐶 1 − 𝑘 21 𝐶 2 (2)
Di  mana 𝑑𝐶 1 𝑑𝑡 adalah  laju  perubahan  konsentrasi  obat dalam sirkulasi sistemik tiap satuan waktu
(mg/L / jam), 𝑘 21 adalah konstanta kecepatan perpindahan obat dari kompartemen pertama (sirkulasi  sistemik)  menuju  ke  kompartemen  kedua atau jaringan  (1/jam)  atau  setara dengan jumlah obat dibagi  dengan  volume  distribusi  pada  kompartemen sentral  ( 𝐷 𝑉 1 ), 𝐶 2 konsentrasi  obat  dalam kompartemen  kedua  (mg/L), 𝑘 12 adalah  konstanta kecepatan perpindahan obat dari kompartemen kedua menuju  ke  kompartemen  pertama, 𝐶 1 adalah konsentrasi  obat  dalam  kompartemen  pertama (mg/L), 𝑘 10 adalah  konstanta  kecepatan  eliminasi obat  (1/jam)  atau  sama  dengan  hasil  kali  antara bersihan  ginjal  dengan  volume  distribusi  pada kompartemen sentral ( 𝐶𝐿 × 𝑉 1 ).
Dengan  cara  yang  sama,  dapat  diturunkan persamaan  diferensial  untuk  model  farmakokinetika dua kompartemen ekstravaskular. Perbedaannya yaitu di  dalam  model  ekstravaskular,  perlu  ditambahkan
parameter-parameter  dalam  proses  absorpsi  dari tempat pemberian menuju ke sirkulasi sistemik. Oleh karena  itu,  persamaan  diferensial  yang  disusun menjadi sedikit berbeda yaitu:
𝑑𝐷 𝑇𝐴 𝑑𝑡 =   −𝑘 𝑎 𝐷 𝑇𝐴 (3) 𝑑𝐶 1 𝑑𝑡 =   −𝑘 10 𝐶 1 − 𝑘 12 𝐶 1 + 𝑘 21 𝐶 2 + 𝐹 𝑟 𝑘 𝑎 𝐷 𝑇𝐴 𝑉 1 (4) 𝑑𝐶 2 𝑑𝑡 =   −𝑘 21 𝐶 2 + 𝑘 12 𝐶 1 (5)
Di mana 𝑑𝐷 𝑇𝐴 𝑑𝑡 adalah laju perubahan konsentrasi obat di  dalam  tempat  absorpsi  (mg/jam), 𝑘 𝑎 adalah konstanta  kecepatan  absorpsi  (1/jam), 𝐷 𝑇𝐴 adalah jumlah  obat  yang  tersedia  dalam  tempat  absorpsi pada  waktu 𝑡 (mg), 𝐹 𝑟 adalah  fraksi  obat  yang terabsorpsi  ke  dalam  sirkulasi  sistemik  yang  dapat diperoleh  dari  perbandingan  antara  luas  daerah  di bawah  kurva  profil  farmakokinetika  ekstravaskular dan  intravena  (tidak  berdimensi  atau  dinyatakan dalam  %), 𝑉 1 adalah  volume  distribusi  pada kompartemen  pertama  atau  sirkulasi  sistemik  (L). Konstanta-konstanta yang menggambarkan kecepatan perpindahan  obat  dari  satu  kompartemen  dari kompartemen yang lain  ini disebut dengan  konstanta mikroskopik.
Persamaan  diferensial  (1) ̶ (5)  sangat penting  untuk  memahami  simulasi  dan  pemodelan farmakokinetika.  Dengan  menggunakan  persamaan diferensial  yang  diprogram  dalam  perangkat  lunak dapat  dilakukan  simulasi  profil  farmakokinetika sehingga  peneliti  maupun  praktisi  klinis  dapat merancang  percobaan  untuk  subjek  uji  atau pengaturan  dosis  terapetik  kepada  pasien  dengan paket  deSolve  (Soetaert et  al ,  2010).  Namun pembahasan  terkait  simulasi  farmakokinetik  berada di  luar  ruang  lingkup  tulisan  ini.  Selain  itu,  tersedia pula  paket-paket  program  yang  mampu  digunakan untuk  melakukan  pencocokan  kurva  ( curve  fitting ) hasil  simulasi  dengan  data  percobaan  untuk menentukan  parameter-parameter  farmakokinetika konvensional  berdasarkan  persamaan  diferensial seperti  WinSAAM  atau  Pkfit  (Lee  &  Lee,  2017; Stefanovski et al , 2003).

## PEMODELAN FARMAKOKINETIKA BERBASIS POPULASI
Farmakokinetika populasi adalah studi farmakokinetika pada tingkatan populasi. Pemodelan farmakokinetik  berbasis  populasi  dapat  dilakukan dengan dua cara yaitu melakukan kalkulasi parameter farmakokinetik dengan menggunakan data kombinasi dari  semua  individu,  menghilangkan  perbedaan individual  (pendekatan  naif  terkelompok)  atau dengan  melakukan  kalkulasi  parameter-parameter farmakokinetika  individual  secara  terpisah  kemudian mengkombinasikan parameter-parameter farmakokinetika dari masing-masing individu sehingga diperoleh parameter-parameter
farmakokinetika  rata-rata  yang  dijadikan  sebagai parameter-parameter  farmakokinetika  populasi  atau disebut  pendekatan  dua  tahap  (Mould  &  Upton, 2012).
Dalam  proses  analisis  data  pada  studi farmakokinetika berbasis populasi, digunakan analisis  pemodelan  efek  campuran.  Pemodelan  efek campuran  pada  dasarnya  adalah  analisis  regresi  baik linier  atau  non-linier  yang  melibatkan  efek  tetap  dan acak.  Baik  efek  tetap  maupun  acak  adalah  koefisien atau  parameter  regresi  yang  menggambarkan hubungan  antara  variabel  bebas  dan  tergantung, hanya  saja  efek  tetap  memodelkan  data  dalam  satu populasi  dan  bernilai  tetap  atau  berubah  sebagai fungsi kovariat, sedangkan efek acak merepresentasikan data dalam sub-kelompok sehingga mempunyai nilai yang bervariasi (Comets et al , 2017; Hamilton, 2013). Dalam konteks farmakokinetika,  digunakan  regresi  non-linier  untuk menggambarkan  hubungan  antara  konsentrasi  obat dalam  plasma  (variabel  tergantung)  dengan  waktu setelah  pemberian  obat  (variabel  bebas)  sedangkan yang  dimaksud  dengan  efek  tetap  dan  acak  adalah parameter-parameter  farmakokinetika  itu  sendiri yang  meliputi  namun  tidak  terbatas  pada  salah  satu atau  lebih  dari: 𝑘 10 , 𝑘 𝑎 , 𝑘 12 , 𝑘 21 , 𝐶𝐿 ,  dan 𝑉 1 tergantung pada model kompartemen yang sesuai dan rute  pemberian  obat,  sedangkan  kovariat  adalah variabel-variabel  lain  di  luar  model  yang  dapat berpengaruh  terhadap  koefisien  regresi  seperti  umur, berat  badan,  atau  gena(de  Alwis et  al ,  1998; Drikvandi, 2017).
Model  efek  campuran  non-linier  dalam  studi farmakokinetika  melibatkan  efek  tetap  yang  terkait dengan parameter-parameter populasi serta efek acak yang  berasal  dari  variabilitas  individu  yang  tidak dapat  dijelaskan  (Tornøe et  al ,  2004).  Variabilitas intra-individual  yang  menggambarkan  perbedaan
antara  nilai  terprediksi  dan  teramati  (residual)  pada suatu  individu  dapat  dikalkulasi  menggunakan persamaan (6) sebagai berikut:
𝑦 𝑖𝑗 = 𝑓 ∅ 𝑖 , 𝑥 𝑖𝑗 + 𝜖 𝑖𝑗 (6)
Di mana i = 1, … , N ; j = 1, … , n i ; y ij adalah respon ke- j untuk individu ke- i ; 𝑓 ∙ adalah fungsi non-linier dari  sebuah  vektor  parameter ∅ yang  spesifik  untuk masing-masing  individu  dan  vektor  prediktor x ij ;  N adalah jumlahindividu, dan n i adalah jumlah pengukuran untuk individu ke- i ; 𝜖 𝑖𝑗 adalah kesalahan residual  yang  diasumsikan  bersifat  independen  dan berupa  variabel  acak  yang  terdistribusi  normal dengan rata-rata nol dan varian konstan sebesar σ 2 . Selanjutnya, model yang dihasilkan menghubungkan  parameter-parameter  dari  individu- individu yang berbeda sebagai berikut:
∅ 𝑖 = 𝐴 𝑖𝑗 𝛽 + 𝐵 𝑖𝑗 𝑏 𝑖 (7)
Di  mana 𝐴 𝑖𝑗 dan 𝐵 𝑖𝑗 adalah  desain  matriks  untuk efek  tetap  vektor 𝛽 dan  efek  acak  vektor 𝑏 𝑖 secara berurutan.  Variasi  inter-individu  dimodelkan  dengan efek  acak  vektor 𝑏 𝑖 yang  terdiri  dari  sejumlah k variabel dengan rata-rata nol yang diasumsikan bebas dan  memiliki  sebaran  yang  dianggap  sebagai distribusi  normal  multivariat  dengan  matriks  varian- kovarian 𝜓 .  Kesalahan  residual 𝜖 𝑖𝑗 dan 𝑏 𝑖 diasumsikan  bebas  untuk  semua i dan j .  Parameter dalam  model  efek  campuran  seperti  yang  tertulis dalam  persamaan  (6)  dan  (7)  dihitung  dengan menggunakan metode maximum likelihood (ML) atau dengan  metode restricted  maximum  likelihood (REML) berdasarkan densitas marjinal dari y .
KALKULASI PARAMETER FARMAKOKINETIKA POPULASI PADA MODEL DUA KOMPARTEMEN

## EKSTRAVASKULAR
Guna  memberikan  pemahaman  yang  lebih dalam  tentang  proses  pemodelan  farmakokinetika berbasis  populasi  dengan  persamaan  diferensial, berikut  ini  disajikan  contoh  kalkulasi  parameter farmakokinetika efavirenz dengan menggunakan data sediaan  referensi  pada  uji  bioekivalensi  yang dipublikasikan  oleh  (Ibarra et  al ,  2016).  Parameter farmakokinetika dikalkulasi menggunakan pendekatan  dua  kompartemen  dengan  persamaan diferensial  seperti  pada  persamaan  (3),  (4),  dan  (5)
dan  dikalkulasi  menggunakan  paket  nlmeODE  dan nlme  (Tornøe et  al ,  2004).  Sistematika  proses kalkulasi parameter farmakokinetika efavirenz dengan  model  dua  kompartemen  pada  tulisan  ini dimulai  dari  penyiapan  data,  penulisan  fungsi diferensial, pemodelan, pengembangan model, pemilihan model terbaik, dan validasi model terpilih.
Sebelum  diolah  dengan  paket  nlmeODE  dan nlme,  data  harus  ditabulasikan  dengan  sistematika tertentu  agar  dapat  diproses  (lihat  contoh  pada Gambar  3).  Selanjutnya,  data  yang  telah  disiapkan dalam  program spreadsheet dapat  disimpan  dalam bentuk  csv,  txt,  atau  xlsx  kemudian  di- import ke dalam  R.  Data  yang  masuk  selanjutnya  diubah menjadi  data  kelompok  dengan  perintah  sebagai berikut:
> library(nlme) > library(nlmeODE) >  data.a  <-  groupedData(conc  ~  Time  | Subject,+ data = as.data.frame(efavirenz))
> data.a$Dose[data.a$Time!=0] <- 0 > data.a$Cmt <- rep(1,dim(data.a)[1]) > par(mar=c(4,4,1,1)) >plot(data.a,
ylab='Konsentrasi dalam
plasma  (mg/L)',  xlab  =  "Waktu  setelah pemberian obat (jam)")

## Gambar 3. Contoh format pemasukan data dengan MS Excel
Dengan  demikian,  diperoleh  data  terkelompok  yang sudah  siap  dianalisis  seperti  yang  terlihat  pada Gambar 4.
Gambar 4. Plot data farmakokinetika efavirenz terkelompok
Setelah  data  siap  untuk  diproses,  dilakukan pemodelan  awal  menggunakan  fungsi  diferensial yang  diturunkan  dari  persamaan  (3)  –  (5).  Hanya saja,  parameter 𝐹 𝑟 tidak  diikutsertakan  dalam perhitungan  karena  membutuhkan  data  intravena. Dalam  pemodelan  ini  digunakan  teknik  pencocokan kurva,    sehingga  diperlukan  parameter  awal  untuk melakukan  perhitungan  yang  dituliskan  dengan perintah start=c() seperti  tertulis  di  bawah.
Parameter awal ini sangat menentukan kebaikan hasil prediksi,  semakin  dekat  parameter  awal  dengan parameter yang sesungguhnya maka proses kalkulasi akan  semakin  baik.  Untuk  mendapatkan  parameter awal  yang  tepat  seringkali  perlu  dilakukan  beberapa kali trial and error .
> twoComp <- list(DiffEq=list( +   dTAdt = ~ -ka*TA,
+   dD1dt = ~ ka*TA-k10*D1-k12*D1+k21*D2, +   dD2dt = ~ k12*D1-k21*D2),
+   ObsEq=list(
+     TA = ~ 0,
+     C1 = ~ D1/V1,
+     C2 = ~ 0), +  Parms=c("ka","k10", "k12", "k21", "V1"), +   States=c("TA","D1","D2"),
+   Init=list(0,0,0))
>  a.model <- nlmeODE(twoComp, data.a, LogParms = T) >  a.nlme  <-  nlme(conc  ~  a.model(ka,  k10, k12, k21, V1, Time, Subject),
+  data  =  data.a,  fixed=ka+k10+k12+k21+V1~1, random = pdDiag(ka+k10+k12+k21+V1~1), +  start=c(ka  =  -0.8,  k10  =  -3.0,  k12  =  - 0.5, k21 = -2.5, V1=3.5),
+ control=list(returnObject=TRUE,msVerbose=TR UE),+ verbose=TRUE)
Setelah  pemodelan  awal  berhasil,  dilakukan pemilihan  model  terbaik  yang  dilakukan  dengan mengembangkan  beberapa  model  yang  mungkin yaitu mengurangi jumlah efek acak serta menambahkan  kovariat  seperti  umur,  berat  badan, dan  status  merokok  satu  per  satu  secara  bertahap ( stepwise ).  Dari  beberapa  kemungkinan  tersebut, dikembangkan  beberapa  model  baik  dengan  kovariat maupun  tanpa  kovariat  kemudian  dilakukan  analisis varian  untuk  mengetahui  model  yang  terbaik.  Model yang  baik  memiliki  nilai Akaike  Information Criterion (AIC)  dan Bayesian  Information  Criterion (BIC)  terendah  serta log-likelihood (LogLik)  yang paling tinggi. Selain itu, model yang sederhana (tidak terlalu  banyak  mengandung  efek  acak  atau  kovariat)
lebih  disukai.  Beberapa  model  alternatf  dapat dikembangkan sebagai berikut:

## #Mengurangi jumlah efek acak
>a1.nlme <- update(a.nlme, random = pdDiag(ka+k12~1)) >a2.nlme <- update(a1.nlme, random = pdDiag(ka~1))
#Apabila  kovariat  berat  badan  dicurigai berpengaruh terhadap V1
> a3.nlme <- update(a1.nlme, fixed = list(ka+k12+k21+k10~1, V1~Wt),  start = c(ka = -0.5, k10 =0.2,  k12 = 1.5, k21 = - 2.5, v1=5, Wt = -0.1))
#Apabila kovariat umur dicurigai
berpengaruh terhadap k10
a4.nlme <- update(a1.nlme, fixed =
list(ka+k12+k21+V1~1, k10~Age),  start = c(ka = 1.0, k10 =-3.0,  k12 = -1.5, k21 = - 2.5, v1=10, Age = -0.01))
Untuk  menguji  kebaikan  fungsi  dari  model-model yang  telah  dikembangkan,  dilakukan  analisis  varian sebagai berikut :
>  anova(a.nlme,  a1.nlme,  a2.nlme,  a3.nlme, a4.nlme)
Berdasarkan  hasil  analisis  varian  dari  model-model yang  telah  dikembangkan  (data  tidak  ditampilkan), dapat  diketahui  bahwa  model  a.nlme  dan  a1.nlme tidak berbeda signifikan. Kedua model ini menghasilkan  nilai  LogLik  paling  tinggi  serta  AIC dan  BIC  paling  rendah,  namun  model  a1.nlme  lebih sederhana  dengan  jumlah  efek  acak  lebih  kecil. Dengan  demikian,  model  a1.nlme  dipilih  sebagai model terbaik untuk menggambarkan farmakokinetika  efavirenz  pada  14  subjek  uji. Selanjutnya,  dilakukan  validasi  terhadap  model terpilih  yang  meliputi  uji  kesesuaian  fungsi  dan  uji asumsi residual.
Kesesuaian fungsi suatu model farmakokinetika  berbasis  populasi  dapat  dinilai berdasarkan  analisis  hubungan  antara  kedekatan nilai-nilai terprediksi dengan data eksperimental yang dapat  diamati  dengan  grafik  prediksi  yang  diperluas dan  plot  hasil  prediksi vs observasi.  Semakin  dekat nilai  terprediksi  dengan  data  eksperimental,  semakin baik kekuatan prediksi dari model yang dikembangkan. Selain itu,  ada beberapa asumsi yang harus terpenuhi, antara lain residual yang ternormalisasi  bersifat  independen  dan  terdistribusi
normal  dengan  rata-rata  nol  dan  varian  sebesar  σ 2 . Selain  itu,  efek  acak  harus  terdistribsusi  normal dengan  rata-rata  nol  dan  matriks  kovarian  sebesar  ψ serta  tidak  tergantung  pada  kesalahan  residual (Tornøe et al , 2004).
Grafik  prediksi  yang  diperluas,  dapat diperoleh  dengan  perintah plot(augPred()) yang akan  menampilkan  plot  antara  data  eksperimental dengan  nilai  prediksi  individu  ( subject )  dan  populasi ( fixed ).  Di  lain  pihak,  untuk  memperoleh  plot konsentrasi  hasil  prediksi  dan  yang  teramati  dapat digunakan perintah plot() . Dalam model farmakokinetika  efavirenz  yang  telah  dikalkulasi, nampak  bahwa  terdapat  kesesuaian  fungsi  yang  baik antara  model  dengan  data  eksperimental  yang ditunjukkan  dengan  kedekatan  antara  data-data ekperimental dengan nilai-nilai terprediksi (Gambar 5 dan Gambar 6).
>plot(augPred(a1.nlme,level=0:1), xlab="Waktu setelah pemberian (jam)", ylab="Konsentrasi dalam plasma (mg/L)") >  plot(a1.nlme,  conc~fitted(.,0),  abline  = c(0,1), xlab='Konsentrasi terprediksi', ylab='Konsentrasi
teramati',
main='Populasi')
>  plot(a1.nlme,  conc~fitted(.,1),  abline  = c(0,1), xlab='Konsentrasi terprediksi', ylab='Konsentrasi teramati', main = "individu")
Selanjutnya,  dalam  uji  diagnostik  residual dapat  diketahui  bahwa  residual  yang  ternormalisasi bersifat  independen  yang  ditunjukkan  dengan  pola sebaran  yang  tidak  membentuk  suatu  pola  tertentu (Gambar  7)  yang  dikuatkan  dengan  plot  ACF  pada Gambar 8 C yang menunjukkan bahwa tidak terdapat autokorelasi  yang  signifikan  pada  taraf  kepercayaan 95% (setelah lag 2, semua nilai autokorelasi berada di bawah  taraf  kepercayaan  yang  digambarkan  dengan garis  putus-putus).  Selain  itu,  residual  yang terstandardisasi mengikuti pola distribusi normal dan memiliki  varian  tetap  yang  ditunjukkan  dengan  plot kuantil  normal  terstandard vs residual  terstandard berbentuk “S” (Gambar 8 A) dan terdistribusi secara simetris  di  sekitar  nol  serta  mempunyai  varian  tetap (Gambar  8C).  Di  lain  pihak,  plot  kuantil  normal terstandard vs efek  acak  menunjukkan  pola  yang menyerupai  normal  (Gambar  9).  Dengan  demikian, dapat disimpulkan bahwa tidak ditemukan pelanggaran  asumsi  yang  cukup  berarti  pada  model yang dipilih.
Nilai  parameter  efek  tetap  dan  efek  acak dapat diperoleh dengan perintah
summary(a1.nlme) dan coef(a1.nlme,augFram
e=T) secara  berurutan.  Meskipun  demikian,  perlu diperhatikan  bahwa  koefisien  yang  diperoleh  berada dalam skala logarima natural. Untuk mengubah skala logaritma  natural  menjadi  bilangan  numerus,  dapat dilakukan transformasi eksponensial dengan  perintah exp() .  Secara  ringkas,  parameter  farmakokinetika efavirenz  hasil  kalkulasi  dapat  ditabulasikan  pada Tabel 1.
Gambar 5. Plot kesesuaian fungsi antara parameter populasi dan individu dengan data eksperimental
Tabel 1. Hasil Perhitungan Parameter
Farmakokinetika  Efavirenz  dengan  paket  nlme  dan nlmeODE Subjek 𝑘 𝑎 𝑘 10 𝑘 12 𝑘 21 𝑉 1 (L) 1 0,169 0,034 0,592 0,051 56027,7 2 0,386 0,034 0,642 0,051 56027,7 3 0,231 0,034 0,596 0,051 56027,7 4 0,524 0,034 0,647 0,051 56027,7 5 0,436 0,034 0,677 0,051 56027,7 6 0,793 0,034 0,547 0,051 56027,7 7 0,665 0,034 0,649 0,051 56027,7 8 0,756 0,034 0,767 0,051 56027,7 9 0,247 0,034 0,510 0,051 56027,7 10 0,442 0,034 0,466 0,051 56027,7 11 0,899 0,034 0,536 0,051 56027,7 12 0,312 0,034 0,318 0,051 56027,7 13 1,021 0,034 0,608 0,051 56027,7 14 0,343 0,034 0,384 0,051 56027,7 Keterangan: 𝑘 10 (1/jam), 𝑘 21 (1/jam), 𝑉 1 (L)  merupakan efek tetap sehingga bernilai konstan untuk seluruh individu sedangkan 𝑘 𝑎 (1/jam) dan 𝑘 12 (1/jam) merupakan efek acak yang bervariasi antar individu.
Apabila dibandingkan dengan paket saemix(Comets et  al ,  2017)  yang  menggunakan algoritma SAEM ( Stochastic Approximation
Expectation  Maximization )  seperti  yang  digunakan dalam Monolix (Chan et al , 2011), teknik pemodelan dengan  persamaan  diferensial  pada  paket  nlme  dan nlmeODE  mempunyai  keunggulan  karena  tidak diperlukan penyelesaian persamaan diferensial dengan  kalkulus  integral  oleh  pengguna.  Hal  ini membuat pengguna mampu mengembangkan pemodelan-pemodelan yang cukup rumit dengan cara yang lebih mudah. Keunggulan lain adalah, nlme dan nlmeODE  mampu  mengkalkulasi  semua  konstanta mikroskopik  secara  langsung  sedangkan  dalam saemix  beberapa  konstanta  mikroskopik  seperti 𝑘 10 dan 𝑘 12 tidak  dapat  dikalkulasi  secara  langsung, tetapi  diturunkan  dari  konstanta  makro 𝜆 1 dan 𝜆 1 (Wijnand,  1988).  Meskipun  demikian,  salah  satu kelemahan  nlme  dan  nlmeODE  dibanding  saemix yaitu  semakin  banyak  parameter  yang  dikalkulasi, semakin banyak pula jumlah data yang dipersyaratkan. Selain  itu, paket nlme dan  nlmeODE ini  tidak  mampu  melakukan  analisis  multi-kovariat sebaik  saemix  dan    tidak  dapat  menampilkan  grafik visual predictive check .
Gambar 6. Diagram pencar antara konsentrasi hasil prediksi dan konsentrasi yang teramati (mg/L) pada parameter farmakokinetika populasi dan individu
> plot(a1.nlme, resid(.,0,type = 'n')~fitted(.,1),id  =  0.05,  abline  =  0, xlab="Konsentrasi terprediksi", ylab="Residual terstandard", main='Populasi') > plot(a1.nlme, resid(.,1,type = 'n')~fitted(.,1),id  =  0.05,  abline  =  0, xlab="Konsentrasi terprediksi", ylab="Residual terstandard", main = "Individu")
Gambar 7. Diagram pencar antara konsentrasi terprediksi (mg/L) dengan residual terstandardisasi pada pemodelan
farmakokinetika populasi dan individu
>  qqnorm(a1.nlme,  ~resid(.,1,type  =  'n'), id  =  0.05,  xlab='Residual  terstandard', ylab='Kuantil normal standard') >  plot(a1.nlme,  Subject~resid(.,1,type  = 'n'), abline = 0, xlab='Residual terstandard', ylab='Subjek') > plot(ACF(a1.nlme), alpha = 0.05, ylab='Autokorelasi')

## Gambar 8. Uji asumsi normalitas dan autokorelasi residual
>  qqnorm(a1.nlme,  ~ranef(.),  xlab=  "Efek acak", ylab='Kuantil normal standard')
Gambar 9. Uji asumsi normalitas efek acak

## KESIMPULAN
Pemodelan farmakokinetika berbasis populasi untuk rute ekstravaskuler dengan model dua kompartemen dapat dilakukan dengan baik menggunakan paket nlme dan nlmeODE yang dibuat berdasarkan  pada  persamaan  diferensial.  Dengan demikian,  paket  nlme  dan  nlmeODE  dapat  dijadikan sebagai  salah  satu  alternatif  perangkat  lunak  sumber terbuka dalam pengembangan model farmakokinetika dua kompartemen berbasis populasi.

## DAFTAR PUSTAKA
Boxenbaum, H. G., Riegelman, S., & Elashoff, R. M. (1974). Statistical estimations in pharmacokinetics. Journal of Pharmacokinetics and Biopharmaceutics , 2 (2), 123–148.
Chan, P. L.  S., Jacqmin, P., Lavielle, M., Mcfadyen,
L., Weatherley, B.(2011). The use of the SAEM algorithm in MONOLIX software for estimation  of  population  pharmacokinetic- pharmacodynamic-viral dynamics parameters of maraviroc  in  asymptomatic  HIV  subjects. J
Pharmacokinet Pharmacodyn , 38 , 41–61. Comets,  E.,  Lavenu,  A.,  &  Lavielle,  M.  (2017). Parameter  Estimation  in  Nonlinear  Mixed Effect Models Using saemix, an R Implementation  of  the  SAEM  Algorithm.
Journal of Statistical Software , 80 (3), 1–41. de  Alwis,  D.  P.,  Aarons,  L.,  &  Palmer,  J.  L.  (1998).
Population  pharmacokinetics  of  ondansetron:  a covariate  analysis. British  Journal  of  Clinical
Pharmacology , 46 (2), 117–125.
Drikvandi,  R.  (2017).  Nonlinear  mixed-effects models  for  pharmacokinetic  data  analysis: assessment  of  the  random-effects  distribution.
Journal of Pharmacokinetics
and Pharmacodynamics , 44 (3), 223–232.
FDA.  (1999).  Guidance  for  Industry  Guidance  for
Industry  Population  Pharmacokinetics,  1–35. Retrieved from http://www.fda.gov/cder/guidance/index.htm Hamilton,  L.  C.  (2013). Statistics  with  STATA (12th ed.). Boston: Brooks/Cole. Ibarra, M., Magallanes,  L., Lorier, M., Vázquez, M., & Fagiolino, P. (2016). Complete dataset for 2- treatment, 2-sequence, 2-period efavirenz bioequivalence  study  conducted  with  nightly dosing. Data in Brief , 7 , 751–754. Jambhekar,  S.  S.,  &  Breen,  P.  J.  (2012). Basic Pharmacokinetics (2nd ed.). Philadelphia. Kovalchik,  S.  (2017).  RISmed:  Download  Content from NCBI Databases. Retrieved from
https://cran.r-project.org/package=RISmed Lee,  C.,  &  Lee,  Y.  (2017).  PKfit :  Data  Analysis Tool for Pharmacokinetics. Mould, D. R., & Upton, R. N. (2012). Basic Concepts in  Population  Modeling,  Simulation,  and Model-Based Drug Development. CPT:
Pharmacometrics  &  Systems  Pharmacology , 1 (9), 1–14.
Pinheiro,  J.,  &  Bates,  D.  M.  (2000). Mixed-Effects
Models in S and S-PLUS . New York: Springer-
Verlag.
R  Core  Team.  (2015).  R:  A  Language  and Environment  for  Statistical  Computing.  R Foundation  for  Statistical  Computing,  Vienna,
Austria. Viena, Austria.
Shargel,  L.,  Yu,  A.  B.  C.,  &  Wu-Pong,  S.  (2012).
Applied biopharmaceutics
and
pharmacokinetics (6th ed.). McGraw-Hill Medical. Retrieved from https://books.google.co.id/books/about/Applied _Biopharmaceutics_Pharmacokinetic.html?id=d
BLreKGy5wkC&redir_esc=y Soetaert,  K.,  Petzoldt,  T.,  &  Setzer,  R.  W.  (2010). Solving  Differential  Equations  in  R :  Package deSolve. Journal  of  Statistical  Software , 33 (9),
1–25.
Stefanovski, D., Moate, P. J., Boston, R. C., & Al., E.
(2003). WinSAAM: a windows-based compartmental  modeling  system. Metabolism:
Clinical and Experimental , 52 (9), 1153–1166.
Tornøe,  C.  W.,  Agersø,  H.,  Jonsson,  E.  N., Madsen,
H., & Nielsen, H. A. (2004). Non-linear mixed- effects pharmacokinetic/pharmacodynamic modelling in NLME using differential equations. Computer Methods and Programs in Biomedicine , 76 (1), 31–40.
Wagner,  J.  G.  (1975).  Application  of  the  Loo-
Riegelman  Absorption  Method. Journal  of
Pharmacokinetics  and  Biopharmaceutics , 3 (1), 51–67. Retrieved from http://hdl.handle.net/2027.42/45059 Wijnand,  H.  P.  (1988).  Pharmacokinetic  model equations  for  the  one-  and  two-compartment models  with  first-order  processes  in  which  the absorption  and  exponential  elimination  or distribution  rate  constants  are  equal. Journal  of
Pharmacokinetics and Biopharmaceutics ,
16 (1), 109–128.
Zhang,  Y.,  Huo,  M.,  Zhou,  J.,  &  Xie,  S.  (2010).
PKSolver: An add-in program for pharmacokinetic  and  pharmacodynamic  data analysis in Microsoft Excel. Computer Methods
and Programs in Biomedicine , 99 (3), 306–314.