diff --git a/R/BD_create.R b/R/BD_create.R
new file mode 100644
index 000000000..6cf2a5d35
--- /dev/null
+++ b/R/BD_create.R
@@ -0,0 +1,415 @@
+# goal is to export a PM model from a file, and return a list of the model parameters
+
+# This will be the example for the workshop
+# 1. create a function to write and save the model to a .json file
+# 2. check the model is passing all the checks in BestDose
+# 3. check the model is running in BestDose (posterior and future predictions)
+
+# ==================================================================____
+# Helper function to perform ODE conversion from model_lib
+# ==================================================================____
+
+#' @title Get the ODE from the model library
+#' @description
+#' `r lifecycle::badge("experimental")`
+#' 
+#' @param model_name A character string for the name of the model in the model library
+#' @return A list of the ODE equations for the model
+
+getODEfromLib <- function(model_name) {
+  # get model library
+  model_lib <- model_lib(show = FALSE)
+  # extract the ODE from the model library
+  ode <- model_lib |>
+    dplyr::filter(Name == model_name) |>
+    dplyr::pull(ode) |>
+    tolower()
+  
+  # get the ODE as a character vector, split by line breaks
+  vect_ode <- unlist(stringr::str_split(ode, "\n"))
+
+  # get Name
+  ode_name <- stringr::str_extract(vect_ode, "(?i)dx\\[[0-9]+\\]")
+  ode_name <- tolower(stringr::str_replace_all(ode_name, "\\[|\\]", ""))
+
+  ode_list <- list()
+  # convert the ODE to a character string
+  for (i in seq_along(vect_ode)) {
+    ode_list[[ode_name[i]]] <- vect_ode[i]
+  }
+  
+  return(ode_list)
+}
+
+
+
+# ==================================================================____
+# Function to extract the model into a BestDose file
+# ==================================================================____
+#' Extract primary from a PM model and format them for BestDose
+#' @param PMmodel A PM model object
+#' @return A list of priors formatted for BestDose
+
+extractPMPrimary <- function(PMmodel) {
+  # extract the priors from the PM model
+  priors <- PMmodel$pri
+  
+  # create a list to store the priors in the format required by BestDose
+  priors_list <- list()
+  
+  # loop through each parameter and extract the min, max, and mean values
+  for (param in names(priors)) {
+    priors_list[[param]] <- list(
+      type = "ab",
+      min = priors[[param]]$min,
+      max = priors[[param]]$max
+    )
+  }
+  
+  return(priors_list)
+}
+
+
+
+#' Extract covariates from a PM model and format them for BestDose
+#' @param PMmodel A PM model object
+#' 
+#' @return A list of covariates formatted for BestDose
+#' 
+
+
+extractPMCovariates <- function(PMmodel) {
+  # early return if covariates are not present in the PM model
+  if (is.null(PMmodel$cov)) {
+    return(NULL)
+  }
+
+  cov <- PMmodel$cov
+
+  # cov names
+  cov_names <- names(cov)
+  
+  cov_list <- list()
+  for (i in seq_along(cov)) {
+    cov_list[[cov_names[i]]] <- list(
+      interp = ifelse(cov[[i]] == 1, "linear", "none")
+    )
+  }
+
+  return(cov_list)
+}
+
+
+#' @title Extract the secondary parameters from a PM model and format them for BestDose
+#' @param PMmodel A PM model object
+#' @return A list of secondary parameters formatted for BestDose
+
+
+extractPMSecondary <- function(PMmodel) {
+  if (is.null(PMmodel$sec)) {
+    return(NULL)
+  }
+
+  sec <- deparse(PMmodel$sec)
+
+  # detect all lines with an assignment operator (= or <-) and trim the white space at the beginning and end of the string
+  bloc <- sec[stringr::str_detect(sec, "=|<-")]
+  bloc <- stringr::str_trim(bloc)
+
+  # get name for the list
+  sec_name <- stringr::str_extract(bloc, "^[a-zA-Z0-9_]+")
+
+  sec_list <- list()
+  for (i in seq_along(bloc)) {
+    sec_list[[sec_name[i]]] <- bloc[i]
+  }
+
+  return(sec_list)
+}
+
+
+
+#' Extract the equations from a PM model and format them for BestDose
+#' @param PMmodel A PM model object
+#' @return A list of equations formatted for BestDose
+
+
+extractPMequation <- function(PMmodel) {
+
+  # equation deparse
+  eqn <- stringr::str_trim(deparse(PMmodel$eqn))
+
+  # get library name
+  library_name <- model_lib(show = FALSE)$Name
+
+  # check if equation are from the library
+  if (any(eqn %in% library_name)) {
+    # extract the model name from the equation and remove NA values
+    model_name <- stringr::str_extract(eqn, paste0("(?i)(", paste(library_name, collapse = "|"), ")"))
+    model_name <- model_name[!is.na(model_name)]
+    return(getODEfromLib(model_name))
+  }
+
+  
+  # capture the derivates equation and trim the white space at the beginning and end of the string
+  # note (?i) is for case insensitive search, and we are looking for "dx[1]", "dx[2]", etc. and capture the whole line
+  bloc <- eqn[stringr::str_detect(eqn, "(?i)dx\\[[0-9]+\\]")]
+
+  # get name for the list
+  eqn_name <- stringr::str_extract(bloc, "dx\\[[0-9]+\\]")
+  eqn_name <- stringr::str_replace_all(eqn_name, "\\[|\\]", "")
+
+  eqn_bloc <- list()
+
+  for (i in seq_along(bloc)) {
+    eqn_bloc[[eqn_name[i]]] <- bloc[i]
+  }
+
+  return(eqn_bloc)
+}
+#' @title Extract the lag time parameters from a PM model and format them for BestDose
+#' @description 
+#' `r lifecycle::badge("experimental")`
+#' 
+#' @param PMmodel A PM model object
+#' @return A list of the lag time parameter if present in the PM model, otherwise NULL
+ 
+
+extractPMLag <- function(PMmodel) {
+  if (is.null(PMmodel$lag)) return(NULL)
+  
+  # extract lag time parameter and format it for BestDose
+  lag <- deparse(PMmodel$lag)
+
+  # extract lag time parameter and trim the white space at the beginning and end of the string
+  # case incensitive search, hence the '(?i)' for "lag[1]", "lag[2]", etc. and capture the whole line
+  bloc <- lag[stringr::str_detect(lag, "(?i)lag\\[[0-9]+\\]")]
+  bloc <- stringr::str_trim(bloc)
+
+  # get name for the list
+  lag_name <- stringr::str_extract(bloc, "(?i)lag\\[[0-9]+\\]")
+  lag_name <- stringr::str_replace_all(lag_name, "\\[|\\]", "")
+
+  lag_list <- list()
+  for (i in seq_along(bloc)) {
+    lag_list[[lag_name[i]]] <- bloc[i]
+  }
+
+  return(lag_list)
+}
+
+
+#' @title Extract the fa parameters from a PM model and format them for BestDose
+#' @description 
+#' `r lifecycle::badge("experimental")`
+#' 
+#' @param PMmodel A PM model object
+#' @return A list of the fa parameter if present in the PM model, otherwise NULL
+
+extractPMFa <- function(PMmodel) {
+if (is.null(PMmodel$fa)) return(NULL)
+
+
+# extract fa parameter and format it for BestDose
+  fa <- deparse(PMmodel$fa)
+
+  # extract fa parameter and trim the white space at the beginning and end of the string
+  # case incensitive search, hence the '(?i)' for "fa[1]", "fa[2]", etc. and capture the whole line
+  bloc <- fa[stringr::str_detect(fa, "(?i)fa\\[[0-9]+\\]")]
+  bloc <- stringr::str_trim(bloc)
+
+  # get name for the list
+  fa_name <- stringr::str_extract(bloc, "(?i)fa\\[[0-9]+\\]")
+  fa_name <- stringr::str_replace_all(fa_name, "\\[|\\]", "")
+
+  fa_list <- list()
+  for (i in seq_along(bloc)) {
+    fa_list[[fa_name[i]]] <- bloc[i]
+  }
+
+  return(fa_list)
+
+}
+
+#' @title Extract the initial parameters from a PM model and format them for BestDose
+#' @description 
+#' `r lifecycle::badge("experimental")`
+#' 
+#' @param PMmodel A PM model object
+#' @return A list of the initial parameters if present in the PM model, otherwise NULL
+
+extractPMInitialVal <- function(PMmodel) {
+if (is.null(PMmodel$ini)) return(NULL)
+
+
+ini <- deparse(PMmodel$ini)
+
+  # extract initial value parameter and trim the white space at the beginning and end of the string
+  # case insensitive search, hence the '(?i)' for "x[1]", "x[2]", etc. and capture the whole line
+  bloc <- ini[stringr::str_detect(ini, "(?i)x\\[[0-9]+\\]")]
+  bloc <- stringr::str_trim(bloc)
+
+  # get name for the list
+  ini_name <- stringr::str_extract(bloc, "(?i)x\\[[0-9]+\\]")
+  ini_name <- stringr::str_replace_all(ini_name, "\\[|\\]", "")
+
+  ini_list <- list()
+  for (i in seq_along(bloc)) {
+    ini_list[[ini_name[i]]] <- bloc[i]
+  }
+
+  return(ini_list)
+}
+
+
+#' @title Extract the output parameters from a PM model and format them for BestDose
+#' @description 
+#' `r lifecycle::badge("experimental")`
+#' 
+#' @param PMmodel A PM model object
+#' @return A list of the output parameter if present in the PM model, otherwise NULL
+
+
+extractPMOuteq <- function(PMmodel) {
+  if (is.null(PMmodel$out)) {
+    cli::cli_abort(c("x" = "The output block is NULL. Please check your PM model."))
+  }
+
+  # deparse the function call
+  out <- deparse(PMmodel$out)
+
+  # extract out time parameter and trim the white space at the beginning and end of the string
+  # case insensitive search, hence the '(?i)' for "y[1]", "y[2]", etc. and capture the whole line
+  bloc <- out[stringr::str_detect(out, "(?i)y\\[[0-9]+\\]")]
+  bloc <- stringr::str_trim(bloc)
+
+  # get name for the list
+  out_name <- stringr::str_extract(bloc, "(?i)y\\[[0-9]+\\]")
+  out_name <- stringr::str_replace_all(out_name, "\\[|\\]", "")
+
+  out_list <- list()
+  for (i in seq_along(bloc)) {
+    out_list[[out_name[i]]] <- bloc[i]
+  }
+
+  return(out_list)
+}
+
+
+#' @title Extract the error parameters from a PM model and format them for BestDose
+#' @description 
+#' `r lifecycle::badge("experimental")`
+#' 
+#' @param PMmodel A PM model object
+#' @return A list of the error parameter if present in the PM model, otherwise NULL
+
+extractPMError <- function(PMmodel) {
+
+  if (is.null(PMmodel$err)) {
+    cli::cli_abort(c("x" = "The error block is NULL. Please check your PM model."))
+  }
+
+  # extract error parameter and format it for BestDose
+  err <- PMmodel$err[[1]]
+  error_list <- list(
+    type = err$type,
+    initial_value = err$initial,
+    coefficient = data.frame(c0 = err$coeff[1], c1 = err$coeff[2], c2 = err$coeff[3], c3 = err$coeff[4])
+  )
+
+  return(error_list)
+}
+
+#' @title Extract the output equations from a PM model and format them for BestDose
+#' @description
+#' `r lifecycle::badge("experimental")`
+#'
+#' @param PM_result A PM model object
+#' @param covariates the covariates information extracted with `extractPMCovariates` function.
+#' @param ode the ordinary differential equations information extracted with `extractPMequation` function.
+#' @param drug_name A character string for the drug name to be used in the description
+#' 
+#' @return A description list for the bestdose model.
+
+createBDdescription <- function(PM_result, covariates = NULL, ode = NULL, drug_name = "Drug") {
+  # gather basic elements of the model to be used in BestDose
+  PMmodel <- PM_result$model$arg_list
+  covariates_list <- NULL
+  # retrieve covariates informations based on the covariates block in the PM model, if it exists
+  if (!is.null(covariates)) {
+    covariates_list <- list(
+      number = length(covariates),
+      names = names(covariates),
+      label = names(covariates),
+      units = rep("to update", length(covariates)),
+      types = rep("numeric", length(covariates)),
+      value = lapply(covariates, function(x) data.frame(min = 0, max = 200, default = 70)),
+      description = lapply(covariates, function(x) "Pls update the description as needed")
+    )
+  }
+
+  # get the number of compartments from the number of ode
+  if (!is.null(ode)) {
+    compartments <- length(ode)
+  } else {
+    compartments <- 1
+  }
+
+  # create the description block
+  description <- list(
+    drug = drug_name,
+    route = data.frame(route = "IV", compartment = 1),
+    name = paste0(drug_name, ".json"),
+    pmx_path = NULL,
+    version = 1.0,
+    compartments = compartments,
+    target = "concentration",
+    target_unit = "mg/L",
+    dose_unit = "mg",
+    description = "This is a model imported from a PM model object.",
+    reference = "This model was imported from a PM model object.",
+    covariates = covariates_list
+  )
+
+  return(description)
+}
+
+#' @title Create a BestDose model list from a PM model object
+#' @description
+#' `r lifecycle::badge("experimental")`
+#' 
+#' @param PM_result A PM model result object
+#' @param drug_name A character string for the drug name to be used in the description
+#' @return A list containing the description, model, and support points formatted for BestDose
+
+createBDmodel <- function(PM_result, drug_name = "Drug") {
+  # gather basic elements of the model to be used in BestDose
+  PMmodel <- PM_result$model$arg_list
+  support_points <- PM_result$final$popPoints
+  ode = extractPMequation(PMmodel)
+
+  # create the model part of the file
+  model_list <- list(
+    primary = extractPMPrimary(PMmodel),
+    covariates = extractPMCovariates(PMmodel),
+    secondary = extractPMSecondary(PMmodel),
+    initial_conditions = extractPMInitialVal(PMmodel),
+    fa = extractPMFa(PMmodel),
+    lag = extractPMLag(PMmodel),
+    equation = ode,
+    out = extractPMOuteq(PMmodel),
+    error = extractPMError(PMmodel)
+  )
+
+  # create the description part of the file
+  description <- createBDdescription(PM_result, covariates = model_list$covariates, ode = ode, drug_name = drug_name)
+
+  # create the R object to write the model file
+  model_file <- list(
+    description = description,
+    model = model_list,
+    support_point = support_points
+  )
+
+  return(model_file)
+}
diff --git a/R/PM_result.R b/R/PM_result.R
index 38b2d5f28..0a45b2602 100755
--- a/R/PM_result.R
+++ b/R/PM_result.R
@@ -301,6 +301,39 @@ PM_result <- R6::R6Class(
     #' @param ... Additional arguments passed to the underlying `fit` method of the model. See the documentation for `fit` for details on supported parameters.
     continue = function(...) {
       self$model$fit(data = self$data, prior = self$final$popPoints, ...)
+    },
+
+    #' @description
+    #' `r lifecycle::badge("experimental")`
+    #' export the model in a JSON file format that can be used in the BestDose Shiny app. 
+    #' This is an experimental feature and may change in future versions.
+    #' 
+    #' @details
+    #' This feature is useful if the user wants to use this specific model to perform MIPD.
+    #' The file will contain most of the features and covariates, but the user will have to specify
+    #' some additional information after importing this file.
+    #' 
+    #' @param drug_name Name of the drug to use in the exported JSON file. Default is "Drug".
+    #' @param dir_path Optional path to save the JSON file. If not provided, the file will be saved in the current working directory with the name "{drug_name}.json".
+    BDexport = function(drug_name = "Drug", dir_path = NULL) {
+      # create model using functions in BD_create.R
+      model_list <- createBDmodel(self, drug_name = drug_name)
+
+      # get the file name to save the JSON file
+      file_name <- if (is.null(dir_path)) {
+        paste0(drug_name, ".json")
+      } else {
+        file.path(dir_path, paste0(drug_name, ".json"))
+      }
+      
+      # #export to JSON
+      jsonlite::toJSON(model_list, pretty = TRUE, auto_unbox = TRUE, null = "null") |>
+      writeLines(con = file_name)
+
+      # provide feedback and return file path invisibly
+      message("Model exported to JSON file: ", file_name)
+      invisible(file_name)
+
     }
   ) # end public
 ) # end PM_result
diff --git a/tests/testthat/test-BD_create.R b/tests/testthat/test-BD_create.R
new file mode 100644
index 000000000..7fabce2fc
--- /dev/null
+++ b/tests/testthat/test-BD_create.R
@@ -0,0 +1,342 @@
+# generate data
+library(Pmetrics)
+PM_result <- NPex
+
+PMmodel <- NPex$model$arg_list
+
+# =================================
+# Test Helper
+# ================================
+testthat::test_that("getODEfromLib returns adequate values", {
+  ode <- getODEfromLib("two_comp_bolus")
+
+  expect_type(ode, "list")
+  expect_equal(length(ode), 2)
+  expect_equal(names(ode), c("dx1", "dx2"))
+  expect_equal(ode$dx1, "dx[1] = b[1] - ka * x[1]")
+  expect_equal(ode$dx2, "dx[2] = r[1] + ka * x[1] - ke * x[2]")
+})
+
+
+
+# =================================
+# Test extractPMPrimary
+# ================================
+
+testthat::test_that("extractPMPrimary", {
+  pri <- extractPMPrimary(PMmodel)
+
+  expect_type(pri, "list")
+  expect_equal(length(pri), 4)
+  expect_equal(names(pri), c("ka", "ke", "v", "tlag1"))
+  expect_equal(names(pri$ka), c("type", "min", "max"))
+})
+
+
+# =================================
+# Test extractPMCovariate
+# ================================
+testthat::test_that("extractPMCovariates return adequate values", {
+  cov <- extractPMCovariates(PMmodel)
+
+  expect_type(cov, "list")
+  expect_equal(length(cov), 5)
+  expect_equal(names(cov), c("wt", "africa", "age", "gender", "height"))
+  expect_equal(cov$wt$interp, "linear")
+  expect_equal(cov$gender$interp, "linear")
+})
+
+testthat::test_that("extractPMCovariates handles NULL", {
+  cov <- NULL
+
+  expect_null(extractPMCovariates(cov))
+})
+
+
+# =================================
+# Test extractPMSecondary
+# ================================
+
+testthat::test_that("extractPMSecondary return adequate with ode", {
+  PMmodel$sec <- function() {
+    v = v0 * wt
+    ke = ke0 * age
+  }
+
+  sec <- extractPMSecondary(PMmodel)
+
+  expect_type(sec, "list")
+  expect_equal(length(sec), 2)
+  expect_equal(names(sec), c("v", "ke"))
+  expect_equal(sec$v, "v = v0 * wt")
+  expect_equal(sec$ke, "ke = ke0 * age")
+})
+
+
+testthat::test_that("extractPMSecondary", {
+  cov <- extractPMSecondary(PMmodel)
+
+  expect_null(cov)
+})
+
+
+
+# =================================
+# Test extractPMequation
+# ================================
+
+testthat::test_that("extractPMequation return adequate with ode", {
+  eqn <- extractPMequation(PMmodel)
+
+  expect_type(eqn, "list")
+  expect_equal(length(eqn), 2)
+  expect_equal(names(eqn), c("dx1", "dx2"))
+  expect_equal(eqn$dx1, "dx[1] = b[1] - ka * x[1]")
+  expect_equal(eqn$dx2, "dx[2] = rateiv[1] + ka * x[1] - ke * x[2]")
+})
+
+testthat::test_that("extractPMequation return adequate values with mod_lib()", {
+  PMmodel$eqn <- function() {
+    two_comp_bolus
+  }
+
+  eqn <- extractPMequation(PMmodel)
+
+  expect_type(eqn, "list")
+  expect_equal(length(eqn), 2)
+  expect_equal(names(eqn), c("dx1", "dx2"))
+  expect_equal(eqn$dx1, "dx[1] = b[1] - ka * x[1]")
+  expect_equal(eqn$dx2, "dx[2] = r[1] + ka * x[1] - ke * x[2]")
+})
+
+
+
+# =================================
+# Test extractPMLag
+# ================================
+
+testthat::test_that("extractPMLag return adequate with ode", {
+  lag <- extractPMLag(PMmodel)
+
+  expect_type(lag, "list")
+  expect_equal(length(lag), 1)
+  expect_equal(names(lag), c("lag1"))
+  expect_equal(lag$lag1, "lag[1] = tlag1")
+})
+
+testthat::test_that("extractPMLag : NULL lag block", {
+  PMmodel$lag <- NULL
+  lag <- extractPMLag(PMmodel)
+
+  expect_null(lag)
+})
+
+# =================================
+# Test extractPMFa
+# ================================
+
+testthat::test_that("extractPMFa handle fa function", {
+  PMmodel$fa <- function() {
+    fa[1] = 0.8
+  }
+
+  fa <- extractPMFa(PMmodel)
+
+  expect_type(fa, "list")
+  expect_equal(length(fa), 1)
+  expect_equal(names(fa), c("fa1"))
+  expect_equal(fa$fa1, "fa[1] = 0.8")
+})
+
+
+testthat::test_that("extractPMFa handle NULL", {
+  fa <- extractPMFa(PMmodel)
+
+  expect_null(fa)
+})
+
+# =================================
+# Test extractPMInitialVal
+# ================================
+
+testthat::test_that("extractPMInitialVal handle NULL", {
+  iv <- extractPMInitialVal(PMmodel)
+
+  expect_null(iv)
+})
+
+testthat::test_that("extractPMInitialVal handle NULL", {
+  PMmodel$ini <- function() {
+    X[1] = 10
+    X[2] = 10 * v
+  }
+
+  iv <- extractPMInitialVal(PMmodel)
+
+  expect_equal(length(iv), 2)
+  expect_equal(names(iv), c("X1", "X2"))
+})
+
+
+# =================================
+# Test extractPMOuteq
+# =================================
+
+testthat::test_that("extractPMOuteq return adequate values", {
+  outeq <- extractPMOuteq(PMmodel)
+
+  expect_type(outeq, "list")
+  expect_equal(length(outeq), 1)
+  expect_equal(names(outeq), "y1")
+  expect_equal(outeq$y1, "y[1] = x[2]/v")
+})
+
+
+testthat::test_that("extractPMOuteq handles NULL", {
+  PMmodel$out <- NULL
+  expect_error(
+    extractPMOuteq(PMmodel),
+    "The output block is NULL. Please check your PM model."
+  )
+})
+
+# =================================
+# Test extractPMError
+# =================================
+
+testthat::test_that("extractPMError return adequate values", {
+  err <- extractPMError(PMmodel)
+
+  expect_type(err, "list")  
+  expect_equal(length(err), 3)
+  expect_equal(names(err), c("type", "initial_value", "coefficient"))
+  expect_equal(err$type, "proportional")
+  expect_equal(err$initial_value, 5)
+  expect_equal(err$coefficient, data.frame(c0 = 0.02, c1 = 0.05, c2 = -0.0002, c3 = 0.0))
+})
+
+testthat::test_that("extractPMError handles NULL", {
+  PMmodel$err <- NULL
+  expect_error(
+    extractPMError(PMmodel),
+    "The error block is NULL. Please check your PM model."
+  )
+})
+
+# =================================
+# Test createBDdescription
+# =================================
+
+testthat::test_that("createBDdescription return adequate values if cov is NULL", {
+  description <- createBDdescription(PM_result)
+  expected_desc_names <- c("drug", "route", "name", "pmx_path","version", "compartments", "target", "target_unit",  "dose_unit", "description", "reference", "covariates")
+  expect_type(description, "list")  
+  expect_equal(length(description), 12)
+  expect_equal(names(description), expected_desc_names)
+  expect_equal(description$drug, "Drug")
+  expect_equal(description$route, data.frame(route = "IV", compartment = 1))
+  expect_equal(description$name, "Drug.json")
+  expect_null(description$pmx_path)
+  expect_equal(description$version, 1.0)
+  expect_equal(description$compartments, 1) 
+  expect_equal(description$target, "concentration")
+  expect_equal(description$target_unit, "mg/L")
+  expect_equal(description$dose_unit, "mg")
+  expect_equal(description$description, "This is a model imported from a PM model object.")
+  expect_equal(description$reference, "This model was imported from a PM model object.")
+  expect_null(description$covariates)
+})
+
+testthat::test_that("createBDdescription return adequate values if any cov is present", {
+  ode <- extractPMequation(PMmodel)
+  covariates <- list(
+    wt = list(interp = "linear"),
+    africa = list(interp = "none"),
+    age = list(interp = "linear")
+  )
+  
+  description <- createBDdescription(PM_result, covariates = covariates, ode = ode)
+  
+  expected_desc_names <- c("drug", "route", "name", "pmx_path","version", "compartments", "target", "target_unit",  "dose_unit", "description", "reference", "covariates")
+  expect_type(description, "list")  
+  expect_equal(length(description), 12)
+  expect_equal(names(description), expected_desc_names)
+  expect_equal(description$drug, "Drug")
+  expect_equal(description$route, data.frame(route = "IV", compartment = 1))
+  expect_equal(description$name, "Drug.json")
+  expect_null(description$pmx_path)
+  expect_equal(description$version, 1.0)
+  expect_equal(description$compartments, 2) # should be 2 here but too many manipulation to get the right value
+  expect_equal(description$target, "concentration")
+  expect_equal(description$target_unit, "mg/L")
+  expect_equal(description$dose_unit, "mg")
+  expect_equal(description$description, "This is a model imported from a PM model object.")
+  expect_equal(description$reference, "This model was imported from a PM model object.")
+
+  # covariate checks
+  expect_type(description$covariates, "list")
+  expect_equal(description$covariates$number, 3)
+  expect_equal(description$covariates$names, c("wt", "africa", "age"))
+  expect_equal(description$covariates$label, c("wt", "africa", "age"))
+  expect_equal(description$covariates$units, c("to update", "to update", "to update"))
+  expect_equal(description$covariates$type, c("numeric", "numeric", "numeric"))
+  
+  # cov value
+  expect_type(description$covariates$value, "list")
+  expect_equal(length(description$covariates$value), 3)
+  expect_equal(names(description$covariates$value), c("wt", "africa", "age"))
+  expect_equal(description$covariates$value$wt, data.frame(min = 0, max = 200, default = 70))
+  expect_equal(description$covariates$value$africa, data.frame(min = 0, max = 200, default = 70))
+  expect_equal(description$covariates$value$age, data.frame(min = 0, max = 200, default = 70))
+
+  # cov description
+  expect_type(description$covariates$description, "list")
+  expect_equal(length(description$covariates$description), 3)
+  expect_equal(names(description$covariates$description), c("wt", "africa", "age"))
+  expect_equal(description$covariates$description$wt, "Pls update the description as needed")
+  expect_equal(description$covariates$description$africa, "Pls update the description as needed")
+  expect_equal(description$covariates$description$age, "Pls update the description as needed")
+})
+
+
+
+# =================================
+# Test createBDmodel
+# =================================
+
+testthat::test_that("createBDmodel return adequate values", {
+  ode <- extractPMequation(PMmodel)
+  bd_mod <- createBDmodel(PM_result)
+
+  expect_type(bd_mod, "list")  
+  expect_equal(length(bd_mod), 3)
+  expect_equal(names(bd_mod), c("description", "model", "support_point"))
+  
+  # check description
+  expect_type(bd_mod$description, "list")
+  expect_equal(length(bd_mod$description), 12)
+  expect_equal(names(bd_mod$description), c("drug", "route", "name", "pmx_path","version", "compartments", "target", "target_unit",  "dose_unit", "description", "reference", "covariates"))
+  expect_equal(bd_mod$description$drug, "Drug")
+  expect_equal(bd_mod$description$route, data.frame(route = "IV", compartment = 1))
+  expect_equal(bd_mod$description$name, "Drug.json")
+  expect_null(bd_mod$description$pmx_path)
+  expect_equal(bd_mod$description$version, 1.0)
+  expect_equal(bd_mod$description$compartments, 2) # should be
+  expect_equal(bd_mod$description$target, "concentration")
+  expect_equal(bd_mod$description$target_unit, "mg/L")
+  expect_equal(bd_mod$description$dose_unit, "mg")
+  expect_equal(bd_mod$description$description, "This is a model imported from a PM model object.")
+  expect_equal(bd_mod$description$reference, "This model was imported from a PM model object.")
+  expect_type(bd_mod$description$covariates, "list")
+
+  # check model
+  expect_type(bd_mod$model, "list")
+  expect_equal(length(bd_mod$model), 9)
+  expect_equal(names(bd_mod$model), c("primary", "covariates", "secondary", "initial_conditions", "fa", "lag", "equation", "out", "error"))
+
+  # check support points
+  expect_type(bd_mod$support_point, "list")
+  expect_equal(ncol(bd_mod$support_point), 5)
+  expect_equal(nrow(bd_mod$support_point), 19)
+  expect_equal(names(bd_mod$support_point), c("ka", "ke", "v", "tlag1", "prob"))
+})