diff --git a/R/BD_create.R b/R/BD_create.R
new file mode 100644
index 000000000..f87893caf
--- /dev/null
+++ b/R/BD_create.R
@@ -0,0 +1,693 @@
+# Export a Pmetrics fit as a BestDose model file (data/drug_model/<name>.json)
+
+# ==================================================================____
+# Helpers
+# ==================================================================____
+
+#' @title BestDose block key
+#' @description Mirrors `parse_equation_lines()` in model_schema.rs: strip everything that is not
+#' alphanumeric and lowercase it, so `dx[1]` -> `dx1`, `Y[2]` -> `y2`, `Ke` -> `ke`.
+#'
+#' @param x A character vector of left-hand sides
+#' @return A character vector of BestDose block keys
+
+bdKey <- function(x) tolower(gsub("[^A-Za-z0-9]", "", x))
+
+#' @title Drop NULL entries from a list
+#' @description BestDose deserializes an absent key to its default, but the `{}` jsonlite writes
+#' for a NULL is a hard error on scalar fields such as `compartment$peripheral`.
+#'
+#' @param x A list
+#' @return The list without its NULL elements
+
+dropNulls <- function(x) Filter(Negate(is.null), x)
+
+#' @title Deparse a PM model block, one line per statement
+#' @description `deparse()` wraps at 60 characters, which splits a long equation across several
+#' elements and silently drops its tail once the lines are filtered by pattern. Deparsing each
+#' statement of the block on its own keeps one equation per line whatever its length.
+#'
+#' @param block A PM model block: a function, a braced expression or a character vector
+#' @return A character vector with one trimmed line per statement
+
+deparseLines <- function(block) {
+  if (is.character(block)) {
+    return(stringr::str_trim(block))
+  }
+
+  expr <- if (is.function(block)) body(block) else block
+  statements <- if (is.expression(expr)) {
+    as.list(expr)
+  } else if (is.call(expr) && identical(expr[[1]], as.name("{"))) {
+    as.list(expr)[-1]
+  } else {
+    list(expr)
+  }
+
+  # a single statement longer than 500 characters still wraps, which no PK equation
+  # reaches. Raise the cutoff, or paste the pieces back, if one ever does.
+  vapply(
+    statements,
+    function(statement) {
+      paste(
+        stringr::str_trim(deparse(statement, width.cutoff = 500L)),
+        collapse = " "
+      )
+    },
+    character(1)
+  )
+}
+
+#' @title Canonical BestDose covariate name
+#' @description model_check.rs rejects the aliases listed in its `COVARIATE_NAME_SYNONYMS`, so
+#' rename them here rather than shipping a file that fails validation on import.
+#'
+#' @param x A character vector of covariate names
+#' @return The canonical lowercase names
+
+bdCovariateName <- function(x) {
+  synonyms <- c(
+    wt = "weight",
+    tbw = "weight",
+    gfr = "crcl",
+    egfr = "crcl",
+    ccr = "crcl",
+    clcr = "crcl",
+    dfg = "crcl",
+    edfg = "crcl",
+    creatinine = "creat",
+    scr = "creat",
+    ht = "height",
+    size = "height",
+    imc = "bmi",
+    body_mass_index = "bmi",
+    sc = "bsa",
+    body_surface_area = "bsa",
+    sexe = "sex",
+    gender = "sex"
+  )
+  x <- tolower(x)
+  unname(ifelse(x %in% names(synonyms), synonyms[x], x))
+}
+
+
+#' @title Extract an indexed block from a PM model
+#' @description One parser for every `keyword[N] = ...` block (`ini`, `fa`, `lag`, `y`), replacing
+#' the four copy-pasted extractors this file used to carry. Keys follow the BestDose convention
+#' (`ini1`, `fa1`, `lag2`, `y1`, ...).
+#'
+#' @param block A PM model block (function or expression), may be NULL
+#' @param keyword The block keyword: "ini", "fa", "lag" or "y"
+#' @return A named list of equation lines, or NULL when the block is absent or empty
+
+extractPMBlock <- function(block, keyword) {
+  if (is.null(block)) {
+    return(NULL)
+  }
+
+  # case insensitive, anchored so only the assignment line is captured, not a reference to it
+  pattern <- paste0("(?i)^", keyword, "\\s*\\[[0-9]+\\]")
+  bloc <- deparseLines(block)
+  bloc <- bloc[stringr::str_detect(bloc, pattern)]
+  if (length(bloc) == 0) {
+    return(NULL)
+  }
+
+  block_list <- as.list(bloc)
+  names(block_list) <- bdKey(stringr::str_extract(bloc, pattern))
+
+  # replace "<-" by "="
+  block_list <- lapply(block_list, function(x) {
+    stringr::str_replace(x, "<-", "=")
+  })
+
+  # input route replacte "rateiv[n]" or "R[n]" by "r[n]" and "Bolus[n]" or "b[n]" by "B[n]"
+  # this should only work on administration
+  block_list <- lapply(block_list, function(x) {
+    x <- stringr::str_replace_all(
+      x,
+      "(?i)rateiv\\s*\\[\\s*([0-9]+)\\s*\\]",
+      "r[\\1]"
+    )
+    x <- stringr::str_replace_all(
+      x,
+      "(?i)R\\s*\\[\\s*([0-9]+)\\s*\\]",
+      "r[\\1]"
+    )
+    x <- stringr::str_replace_all(
+      x,
+      "(?i)bolus\\s*\\[\\s*([0-9]+)\\s*\\]",
+      "B[\\1]"
+    )
+    x <- stringr::str_replace_all(
+      x,
+      "(?i)b\\s*\\[\\s*([0-9]+)\\s*\\]",
+      "B[\\1]"
+    )
+    x
+  })
+
+  return(block_list)
+}
+
+# ==================================================================____
+# model block
+# ==================================================================____
+
+#' @title Extract primary parameters from a PM model and format them for BestDose
+#' @param PMmodel A PM model object
+#' @return A named list of `{type, min, max}` priors
+
+extractPMPrimary <- function(PMmodel) {
+  priors <- PMmodel$pri
+  priors_list <- list()
+
+  for (param in names(priors)) {
+    prior <- priors[[param]]
+    # BestDose stores both prior flavours in min/max: "ab" = bounds, "msd" = mean/SD
+    priors_list[[param]] <- if (!is.null(prior$min)) {
+      list(type = "ab", min = prior$min, max = prior$max)
+    } else {
+      list(type = "msd", min = prior$mean, max = prior$sd)
+    }
+  }
+
+  return(priors_list)
+}
+
+
+#' @title Extract covariates from a PM model and format them for BestDose
+#' @param PMmodel A PM model object
+#' @return A named list of `{interp}` settings, or NULL when the model has no covariates
+
+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 <- bdCovariateName(names(cov))
+  if (!identical(cov_names, tolower(names(cov)))) {
+    message(
+      "Covariates renamed to their BestDose canonical name: ",
+      paste(names(cov), "->", cov_names, collapse = ", ")
+    )
+    message(
+      "The secondary/ODE equations still use the original names, update them before saving."
+    )
+  }
+
+  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 named list of secondary equations, or NULL when the block is absent
+
+extractPMSecondary <- function(PMmodel) {
+  if (is.null(PMmodel$sec)) {
+    return(NULL)
+  }
+
+  sec <- deparseLines(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, "=|<-")]
+  # keep the "name = expression" lines only, dropping the function() header deparse adds
+  bloc <- bloc[stringr::str_detect(bloc, "^[a-zA-Z][a-zA-Z0-9_]*\\s*(=|<-)")]
+  if (length(bloc) == 0) {
+    return(NULL)
+  }
+
+  sec_list <- as.list(bloc)
+  names(sec_list) <- bdKey(stringr::str_extract(bloc, "^[a-zA-Z0-9_]+"))
+
+  return(sec_list)
+}
+
+#' @title Extract the error model from a PM model and format it for BestDose
+#' @param PMmodel A PM model object
+#' @return The BestDose error block
+
+extractPMError <- function(PMmodel) {
+  if (is.null(PMmodel$err)) {
+    stop("The error block is NULL. Please check your PM model.")
+  }
+
+  err <- PMmodel$err[[1]]
+  # BestDose only knows "additive" and "proportional" (Pmetrics' L and G)
+  err_type <- switch(
+    tolower(as.character(err$type)),
+    "l" = "additive",
+    "additive" = "additive",
+    "g" = "proportional",
+    "proportional" = "proportional",
+    stop("Unsupported error model type: ", err$type)
+  )
+
+  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)
+}
+
+# ==================================================================____
+# description block
+# ==================================================================____
+
+#' @title Derive the administration routes from the equations
+#' @description BestDose needs at least one `{route, compartment}` entry, and that compartment
+#' index is what the engine actually doses into. Infusion tokens (`R[N]`, `RATEIV[N]`) map to IV,
+#' bolus tokens (`B[N]`) to an extravascular depot. The compartment is inferred, the human label is
+#' a guess: confirm it.
+#'
+#' @param equations The equation list from `extractPMequation`
+#' @return A data.frame of routes
+
+extractBDroutes <- function(equations) {
+  text <- tolower(paste(unlist(equations), collapse = " "))
+  indices <- function(pattern) {
+    unique(as.integer(stringr::str_match_all(text, pattern)[[1]][, 2]))
+  }
+
+  iv <- indices("(?:^|[^a-z0-9_])(?:r|rateiv)\\s*\\[\\s*([0-9]+)\\s*\\]")
+  po <- indices("(?:^|[^a-z0-9_])b\\s*\\[\\s*([0-9]+)\\s*\\]")
+
+  if (length(iv) + length(po) == 0) {
+    message(
+      "No dose input (R[N]/B[N]) found in the equations, defaulting to IV into compartment 1."
+    )
+    return(data.frame(route = "IV", compartment = 1L))
+  }
+
+  routes <- data.frame(
+    route = c(rep("IV", length(iv)), rep("PO", length(po))),
+    compartment = as.integer(c(iv, po))
+  )
+  message(
+    "Routes inferred from the equations, check the labels: ",
+    paste(routes$route, "->", routes$compartment, collapse = ", ")
+  )
+
+  return(routes)
+}
+
+#' @title Build the compartment block
+#' @description `number` counts the *named* compartments (one per output equation), not the ODE
+#' states: an absorption depot has an equation but no output. `outputs[i]$equation` is a Y index.
+#'
+#' @param out The output list from `extractPMBlock(PMmodel$out, "y")`
+#' @return The BestDose compartment block
+
+createBDcompartment <- function(out) {
+  n <- length(out)
+  if (n == 0) {
+    stop("The output block is NULL. Please check your PM model.")
+  }
+
+  # the first output must be named "central", the others follow the BestDose naming convention
+  names_out <- c(
+    "central",
+    "peripheral",
+    "effect",
+    paste0("compartment", seq_len(max(n - 3, 0)) + 3)
+  )[seq_len(n)]
+
+  dropNulls(list(
+    number = n,
+    central = 1L,
+    peripheral = if (n >= 2) 2L else NULL,
+    outputs = data.frame(name = names_out, equation = seq_len(n))
+  ))
+}
+
+# Starting ranges for the covariates BestDose knows about. These are the clinical bounds shown in
+# the UI, not statistical ones: a Pmetrics fit carries none of this, so anything missing here is
+# filled with a placeholder that needs a review.
+BD_COVARIATE_INFO <- list(
+  weight = list(
+    label = "Weight",
+    unit = "kg",
+    min = 2,
+    max = 650,
+    default = 70,
+    description = "Patient weight in kilograms"
+  ),
+  height = list(
+    label = "Height",
+    unit = "cm",
+    min = 30,
+    max = 250,
+    default = 170,
+    description = "Patient height in centimeters"
+  ),
+  age = list(
+    label = "Age",
+    unit = "years",
+    min = 0,
+    max = 120,
+    default = 50,
+    description = "Patient age in years"
+  ),
+  crcl = list(
+    label = "Creatinine clearance",
+    unit = "mL/min",
+    min = 0,
+    max = 200,
+    default = 90,
+    description = "Creatinine clearance in mL/min calculated using the Cockcroft-Gault equation"
+  ),
+  creat = list(
+    label = "Serum creatinine",
+    unit = "mg/dL",
+    min = 0,
+    max = 20,
+    default = 1,
+    description = "Serum creatinine concentration"
+  ),
+  bmi = list(
+    label = "Body mass index",
+    unit = "kg/m2",
+    min = 10,
+    max = 60,
+    default = 25,
+    description = "Body mass index"
+  ),
+  bsa = list(
+    label = "Body surface area",
+    unit = "m2",
+    min = 0.1,
+    max = 3,
+    default = 1.8,
+    description = "Body surface area"
+  )
+)
+
+#' @title Build the covariates description block
+#' @description Only the covariate *names* come from the fit: labels, units, ranges and
+#' descriptions are clinical metadata, taken from `BD_COVARIATE_INFO` when known and filled with a
+#' placeholder otherwise. `I()` keeps the single-covariate vectors as JSON arrays.
+#'
+#' @param covariates The covariate list from `extractPMCovariates`
+#' @return The BestDose covariates description block, or NULL when the model has no covariate
+
+createBDcovariates <- function(covariates) {
+  if (length(covariates) == 0) {
+    return(NULL)
+  }
+
+  cov_names <- names(covariates)
+  info <- lapply(cov_names, function(name) {
+    known <- BD_COVARIATE_INFO[[name]]
+    if (is.null(known)) {
+      message(
+        "Unknown covariate '",
+        name,
+        "': update its label, unit, range and description before using the model."
+      )
+      known <- list(
+        label = name,
+        unit = "to update",
+        min = 0,
+        max = 200,
+        default = 70,
+        description = "Please update the description as needed"
+      )
+    }
+    known
+  })
+  names(info) <- cov_names
+
+  list(
+    number = length(cov_names),
+    names = I(cov_names),
+    label = I(vapply(
+      info,
+      function(x) x$label,
+      character(1),
+      USE.NAMES = FALSE
+    )),
+    units = I(vapply(
+      info,
+      function(x) x$unit,
+      character(1),
+      USE.NAMES = FALSE
+    )),
+    types = I(rep("numeric", length(cov_names))),
+    value = lapply(info, function(x) {
+      data.frame(min = x$min, max = x$max, default = x$default)
+    }),
+    description = lapply(info, function(x) x$description)
+  )
+}
+
+#' @title Build the description block of a BestDose model file
+#'
+#' @param model_list The model block built by `createBDmodel`
+#' @param drug_name The drug the model attaches to (stored lowercase)
+#' @param model_name The model file name, without the .json extension
+#' @param description Free text description of the model
+#' @param reference Citation for the model
+#' @param reference_url URL for the citation
+#'
+#' @return The BestDose description block
+
+createBDdescription <- function(
+  model_list,
+  drug_name,
+  model_name,
+  description = "",
+  reference = "",
+  reference_url = ""
+) {
+  dropNulls(list(
+    drug = tolower(drug_name),
+    route = extractBDroutes(model_list$equation),
+    name = paste0(model_name, ".json"),
+    compartment = createBDcompartment(model_list$out),
+    version = 1L,
+    description = description,
+    reference = reference,
+    reference_url = reference_url,
+    covariates = createBDcovariates(model_list$covariates)
+  ))
+}
+
+# ==================================================================____
+# support points
+# ==================================================================____
+
+#' @title Extract the final cycle population points as BestDose support points
+#' @description One row per support point, one column per primary parameter plus `prob`. BestDose
+#' requires those columns to match the primary parameter names exactly.
+#'
+#' @param PM_result A PM result object
+#' @param primary_names The names of the primary parameters
+#' @return A data.frame of support points
+
+extractBDsupportPoints <- function(PM_result, primary_names) {
+  points <- PM_result$final$popPoints
+  if (is.null(points)) {
+    warning(
+      "No population points found in the PM result, the model file will carry no support point."
+    )
+    return(list())
+  }
+
+  points <- as.data.frame(points)
+  names(points)[tolower(names(points)) %in% c("prob", "probability")] <- "prob"
+
+  missing <- setdiff(primary_names, names(points))
+  if (length(missing) > 0) {
+    stop(
+      "Support points are missing the primary parameter(s): ",
+      paste(missing, collapse = ", ")
+    )
+  }
+  if (!"prob" %in% names(points)) {
+    stop("Support points are missing the 'prob' column.")
+  }
+
+  points[, c(primary_names, "prob"), drop = FALSE]
+}
+
+# ==================================================================____
+# assemble, check and save
+# ==================================================================____
+
+#' @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 The drug the model attaches to
+#' @param model_name The model file name, without the .json extension
+#' @param description Free text description of the model
+#' @param reference Citation for the model
+#' @param reference_url URL for the citation
+#' @return A list containing the description, model and support points formatted for BestDose
+#'
+#' @export
+
+createBDmodel <- function(
+  PM_result,
+  drug_name = "Drug",
+  model_name = drug_name,
+  description = "Model imported from a Pmetrics fit.",
+  reference = "",
+  reference_url = ""
+) {
+  # gather basic elements of the model to be used in BestDose
+  PMmodel <- PM_result$model$arg_list
+
+  # create the model part of the file
+  model_list <- dropNulls(list(
+    primary = extractPMPrimary(PMmodel),
+    covariates = extractPMCovariates(PMmodel),
+    secondary = extractPMSecondary(PMmodel),
+    initial_conditions = extractPMBlock(PMmodel$ini, "ini"),
+    fa = extractPMBlock(PMmodel$fa, "fa"),
+    lag = extractPMBlock(PMmodel$lag, "lag"),
+    equation = extractPMBlock(PMmodel$eqn, "dx"),
+    out = extractPMBlock(PMmodel$out, "y"),
+    error = extractPMError(PMmodel)
+  ))
+
+  model_file <- list(
+    description = createBDdescription(
+      model_list,
+      drug_name,
+      model_name,
+      description,
+      reference,
+      reference_url
+    ),
+    model = model_list,
+    support_point = extractBDsupportPoints(PM_result, names(model_list$primary))
+  )
+
+  checkBDmodel(model_file)
+
+  return(model_file)
+}
+
+#' @title Check a BestDose model list before writing it
+#' @description The blocking subset of `model_check.rs`, i.e. everything that makes BestDose refuse
+#' the file outright. The app runs the full check, warnings included, on import.
+#'
+#' @param model_file A model list built by `createBDmodel`
+#' @return `invisible(TRUE)`, or stops with the list of problems
+
+checkBDmodel <- function(model_file) {
+  description <- model_file$description
+  model <- model_file$model
+  issues <- character(0)
+
+  if (nchar(description$drug) == 0) {
+    issues <- c(issues, "A drug must be assigned to the model.")
+  }
+  if (nchar(description$name) == 0) {
+    issues <- c(issues, "The model name is required.")
+  }
+  if (nrow(description$route) == 0) {
+    issues <- c(issues, "At least one administration route is required.")
+  }
+  if (length(model$primary) < 2) {
+    issues <- c(issues, "At least two primary parameters are required.")
+  }
+
+  # compartment and output consistency
+  n_compartment <- description$compartment$number
+  if (n_compartment < 1 || n_compartment > 6) {
+    issues <- c(issues, "Compartment count must be between 1 and 6.")
+  }
+  if (n_compartment > length(model$equation)) {
+    issues <- c(
+      issues,
+      "Compartment count exceeds the number of differential equations."
+    )
+  }
+  if (length(model$out) != n_compartment) {
+    issues <- c(
+      issues,
+      "The number of output equations does not match the number of named compartments."
+    )
+  }
+  if (
+    !setequal(
+      names(model$out),
+      paste0("y", description$compartment$outputs$equation)
+    )
+  ) {
+    issues <- c(
+      issues,
+      "The output equations do not match the compartment output indices."
+    )
+  }
+
+  # a route may dose into an implicit depot, so it is bounded by the ODE count, not the compartment count
+  state_count <- max(length(model$equation), n_compartment)
+  if (
+    any(
+      description$route$compartment < 1 |
+        description$route$compartment > state_count
+    )
+  ) {
+    issues <- c(
+      issues,
+      "A route maps to a compartment outside the model's range."
+    )
+  }
+
+  # every covariate declared in the description needs a model level interpolation, and the reverse
+  cov_described <- if (is.null(description$covariates)) {
+    character(0)
+  } else {
+    as.character(description$covariates$names)
+  }
+  if (!setequal(cov_described, names(model$covariates))) {
+    issues <- c(
+      issues,
+      "The covariates in the description and in the model block do not match."
+    )
+  }
+
+  # support point columns must be exactly the primary parameters plus prob
+  if (
+    length(model_file$support_point) > 0 &&
+      !setequal(
+        names(model_file$support_point),
+        c(names(model$primary), "prob")
+      )
+  ) {
+    issues <- c(
+      issues,
+      "The support point columns do not match the primary parameters."
+    )
+  }
+
+  if (length(issues) > 0) {
+    stop(
+      "The model is not valid for BestDose:\n- ",
+      paste(issues, collapse = "\n- ")
+    )
+  }
+
+  invisible(TRUE)
+}
diff --git a/R/PM_result.R b/R/PM_result.R
index 03849b0ba..6e6d02ef2 100755
--- a/R/PM_result.R
+++ b/R/PM_result.R
@@ -299,6 +299,36 @@ 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. 
+    #' This is an experimental feature and may change in future versions.
+    #' 
+    #' @details
+    #' This feature is useful if the user want to use this specific model to perform MIPD.
+    #' The file will contains 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".
+    
+    export_bestdose = function(drug_name = "Drug", dir_path = NULL) {
+      # create model using functions in BD_create.R
+      model_file <- createBDmodel(self, drug_name = drug_name)
+      path = ifelse(is.null(dir_path), paste0(drug_name, ".json"), paste0(dir_path, "/", drug_name, ".json"))
+      
+      # #export to JSON
+      jsonlite::write_json(
+        model_file,
+        path,
+        auto_unbox = TRUE,
+        dataframe = "rows",
+        null = "null",
+        pretty = TRUE
+      )
+
     }
   ) # 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..1bccee7d0
--- /dev/null
+++ b/tests/testthat/test-BD_create.R
@@ -0,0 +1,476 @@
+# Checks for BD_create.R — run with testthat after sourcing the exporter:
+#   source("docs/BD_create.R"); testthat::test_file("docs/test-BD_create.R")
+
+# generate data
+library(Pmetrics)
+library(testthat)
+PM_result <- NPex
+
+PMmodel <- NPex$model$arg_list
+
+# =================================
+# Test Helper
+# ================================
+
+test_that("bdKey strips brackets and lowercases", {
+  expect_equal(
+    bdKey(c("dx[1]", "Y[2]", "Ke", "lag[10]")),
+    c("dx1", "y2", "ke", "lag10")
+  )
+})
+
+test_that("dropNulls removes the absent blocks", {
+  expect_equal(dropNulls(list(a = 1, b = NULL, c = "x")), list(a = 1, c = "x"))
+})
+
+test_that("deparseLines keeps one statement per line whatever its length", {
+  block <- function() {
+    dx[1] <- r[1] -
+      ke * x[1] -
+      kcp * x[1] +
+      kpc * x[2] -
+      kcp2 * x[1] +
+      kpc2 * x[3] -
+      kcp3 * x[1]
+    dx[2] <- kcp * x[1] - kpc * x[2]
+  }
+
+  lines <- deparseLines(block)
+
+  expect_length(lines, 2)
+  expect_equal(
+    lines[1],
+    "dx[1] <- r[1] - ke * x[1] - kcp * x[1] + kpc * x[2] - kcp2 * x[1] + kpc2 * x[3] - kcp3 * x[1]"
+  )
+  expect_equal(lines[2], "dx[2] <- kcp * x[1] - kpc * x[2]")
+})
+
+
+test_that("extractPMBlock reads any indexed block and handles NULL", {
+  lag <- extractPMBlock(PMmodel$lag, "lag")
+
+  expect_type(lag, "list")
+  expect_equal(names(lag), "lag1")
+  expect_equal(lag$lag1, "lag[1] = tlag1")
+
+  expect_null(extractPMBlock(NULL, "fa"))
+  expect_null(extractPMBlock(
+    function() {
+      x <- 1
+    },
+    "fa"
+  ))
+})
+
+# =================================
+# Test extractPMPrimary
+# ================================
+
+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"))
+  expect_equal(pri$ka$type, "ab")
+  expect_equal(pri$ka$min, 0.1)
+  expect_equal(pri$ka$max, 0.9)
+})
+
+# =================================
+# Test extractPMCovariate
+# ================================
+
+test_that("bdCovariateName maps the aliases BestDose rejects", {
+  expect_equal(
+    bdCovariateName(c("wt", "GFR", "gender", "age")),
+    c("weight", "crcl", "sex", "age")
+  )
+})
+
+test_that("extractPMCovariates return adequate values", {
+  cov <- extractPMCovariates(PMmodel)
+
+  expect_type(cov, "list")
+  expect_equal(length(cov), 5)
+  # wt and gender are aliases, renamed to the canonical BestDose names
+  expect_equal(names(cov), c("weight", "africa", "age", "sex", "height"))
+  expect_equal(cov$weight$interp, "linear")
+  expect_equal(cov$sex$interp, "linear")
+})
+
+test_that("extractPMCovariates handles NULL", {
+  cov <- NULL
+
+  expect_null(extractPMCovariates(cov))
+})
+
+# =================================
+# Test extractPMSecondary
+# ================================
+
+test_that("extractPMSecondary return adequate with ode", {
+  pm <- PMmodel
+  pm$sec <- function() {
+    v <- v0 * wt
+    ke <- ke0 * age
+  }
+
+  sec <- extractPMSecondary(pm)
+
+  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")
+})
+
+test_that("extractPMSecondary", {
+  cov <- extractPMSecondary(PMmodel)
+
+  expect_null(cov)
+})
+
+
+# =================================
+# Test extractPMBlock on fa / ini / out
+# =================================
+
+test_that("extractPMBlock handle fa function", {
+  pm <- PMmodel
+  pm$fa <- function() {
+    fa[1] <- 0.8
+  }
+
+  fa <- extractPMBlock(pm$fa, "fa")
+
+  expect_type(fa, "list")
+  expect_equal(length(fa), 1)
+  expect_equal(names(fa), c("fa1"))
+  expect_equal(fa$fa1, "fa[1] = 0.8")
+})
+
+test_that("extractPMBlock handle NULL fa and ini", {
+  expect_null(extractPMBlock(PMmodel$fa, "fa"))
+  expect_null(extractPMBlock(PMmodel$ini, "ini"))
+})
+
+test_that("extractPMBlock return adequate output equations", {
+  outeq <- extractPMBlock(PMmodel$out, "y")
+
+  expect_type(outeq, "list")
+  expect_equal(length(outeq), 1)
+  expect_equal(names(outeq), "y1")
+  expect_equal(outeq$y1, "y[1] = x[2]/v")
+})
+
+# =================================
+# Test extractPMError
+# =================================
+
+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.0, c1 = 0.1, c2 = 0.0, c3 = 0.0)
+  )
+})
+
+test_that("extractPMError translates the Pmetrics L/G types", {
+  pm <- PMmodel
+  pm$err <- list(list(type = "L", initial = 1, coeff = c(0.1, 0.2, 0, 0)))
+  expect_equal(extractPMError(pm)$type, "additive")
+
+  pm$err <- list(list(
+    type = "unknown",
+    initial = 1,
+    coeff = c(0.1, 0.2, 0, 0)
+  ))
+  expect_error(extractPMError(pm), "Unsupported error model type")
+})
+
+# =================================
+# Test extractBDroutes
+# =================================
+
+test_that("extractBDroutes falls back to IV in compartment 1", {
+  routes <- extractBDroutes(list(dx1 = "dx[1] = -ke * x[1]"))
+
+  expect_equal(nrow(routes), 1)
+  expect_equal(routes$route, "IV")
+  expect_equal(routes$compartment, 1L)
+})
+
+# =================================
+# Test createBDcompartment
+# =================================
+
+test_that("createBDcompartment names the outputs and drops the absent peripheral", {
+  one <- createBDcompartment(list(y1 = "y[1] = x[1]/v"))
+
+  expect_equal(names(one), c("number", "central", "outputs"))
+  expect_equal(one$number, 1)
+  expect_equal(one$outputs, data.frame(name = "central", equation = 1L))
+
+  two <- createBDcompartment(list(y1 = "y[1] = x[1]/v", y2 = "y[2] = x[2]/v"))
+
+  expect_equal(names(two), c("number", "central", "peripheral", "outputs"))
+  expect_equal(two$peripheral, 2L)
+  expect_equal(two$outputs$name, c("central", "peripheral"))
+  expect_equal(two$outputs$equation, c(1L, 2L))
+})
+
+test_that("createBDcompartment stops without an output block", {
+  expect_error(createBDcompartment(NULL), "output block is NULL")
+})
+
+# =================================
+# Test createBDcovariates
+# =================================
+
+test_that("createBDcovariates fills the known ranges and flags the others", {
+  covariates <- list(
+    weight = list(interp = "linear"),
+    africa = list(interp = "none"),
+    age = list(interp = "linear")
+  )
+
+  cov_desc <- createBDcovariates(covariates)
+
+  expect_type(cov_desc, "list")
+  expect_equal(cov_desc$number, 3)
+  expect_equal(as.character(cov_desc$names), c("weight", "africa", "age"))
+  expect_equal(as.character(cov_desc$label), c("Weight", "africa", "Age"))
+  expect_equal(as.character(cov_desc$units), c("kg", "to update", "years"))
+  expect_equal(as.character(cov_desc$types), c("numeric", "numeric", "numeric"))
+
+  # cov value
+  expect_type(cov_desc$value, "list")
+  expect_equal(names(cov_desc$value), c("weight", "africa", "age"))
+  expect_equal(
+    cov_desc$value$weight,
+    data.frame(min = 2, max = 650, default = 70)
+  )
+  expect_equal(
+    cov_desc$value$africa,
+    data.frame(min = 0, max = 200, default = 70)
+  )
+  expect_equal(cov_desc$value$age, data.frame(min = 0, max = 120, default = 50))
+
+  # cov description
+  expect_equal(names(cov_desc$description), c("weight", "africa", "age"))
+  expect_equal(cov_desc$description$weight, "Patient weight in kilograms")
+  expect_equal(
+    cov_desc$description$africa,
+    "Please update the description as needed"
+  )
+})
+
+test_that("createBDcovariates returns NULL without covariates", {
+  expect_null(createBDcovariates(NULL))
+  expect_null(createBDcovariates(list()))
+})
+
+# =================================
+# Test createBDdescription
+# =================================
+
+test_that("createBDdescription return adequate values if cov is NULL", {
+  model_list <- list(
+    equation = extractPMBlock(PMmodel$eqn, "dx"),
+    out = extractPMBlock(PMmodel$out, "y")
+  )
+
+  description <- createBDdescription(model_list, "Amikacin", "AmikacinPM")
+
+  expected_desc_names <- c(
+    "drug",
+    "route",
+    "name",
+    "compartment",
+    "version",
+    "description",
+    "reference",
+    "reference_url"
+  )
+  expect_type(description, "list")
+  expect_equal(names(description), expected_desc_names)
+  expect_equal(description$drug, "amikacin")
+  expect_equal(
+    description$route,
+    data.frame(route = c("IV", "PO"), compartment = c(1L, 1L))
+  )
+  expect_equal(description$name, "AmikacinPM.json")
+  expect_equal(description$version, 1L)
+  expect_equal(description$compartment$number, 1)
+  expect_equal(description$description, "")
+  expect_equal(description$reference, "")
+  expect_equal(description$reference_url, "")
+  expect_null(description$covariates)
+})
+
+test_that("createBDdescription return adequate values if any cov is present", {
+  model_list <- list(
+    covariates = extractPMCovariates(PMmodel),
+    equation = extractPMBlock(PMmodel$eqn, "dx"),
+    out = extractPMBlock(PMmodel$out, "y")
+  )
+
+  description <- createBDdescription(
+    model_list,
+    "Amikacin",
+    "AmikacinPM",
+    description = "a model",
+    reference = "a citation",
+    reference_url = "https://example.org"
+  )
+
+  expected_desc_names <- c(
+    "drug",
+    "route",
+    "name",
+    "compartment",
+    "version",
+    "description",
+    "reference",
+    "reference_url",
+    "covariates"
+  )
+  expect_equal(names(description), expected_desc_names)
+  expect_equal(description$description, "a model")
+  expect_equal(description$reference, "a citation")
+  expect_equal(description$reference_url, "https://example.org")
+
+  # covariate checks
+  expect_type(description$covariates, "list")
+  expect_equal(description$covariates$number, 5)
+  expect_equal(
+    as.character(description$covariates$names),
+    c("weight", "africa", "age", "sex", "height")
+  )
+})
+
+# =================================
+# Test extractBDsupportPoints
+# =================================
+
+test_that("extractBDsupportPoints keeps the primary parameters and prob", {
+  points <- extractBDsupportPoints(PM_result, c("ka", "ke", "v", "tlag1"))
+
+  expect_s3_class(points, "data.frame")
+  expect_equal(nrow(points), 17)
+  expect_equal(names(points), c("ka", "ke", "v", "tlag1", "prob"))
+})
+
+test_that("extractBDsupportPoints stops on a missing parameter", {
+  expect_error(
+    extractBDsupportPoints(PM_result, c("ka", "cl")),
+    "missing the primary parameter"
+  )
+})
+
+# =================================
+# Test createBDmodel
+# =================================
+
+test_that("createBDmodel return adequate values", {
+  bd_mod <- createBDmodel(
+    PM_result,
+    drug_name = "Amikacin",
+    model_name = "AmikacinPM"
+  )
+
+  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(
+    names(bd_mod$description),
+    c(
+      "drug",
+      "route",
+      "name",
+      "compartment",
+      "version",
+      "description",
+      "reference",
+      "reference_url",
+      "covariates"
+    )
+  )
+  expect_equal(bd_mod$description$drug, "amikacin")
+  expect_equal(bd_mod$description$name, "AmikacinPM.json")
+  expect_equal(bd_mod$description$version, 1L)
+  expect_equal(bd_mod$description$compartment$number, 1)
+  expect_equal(
+    bd_mod$description$compartment$outputs,
+    data.frame(name = "central", equation = 1L)
+  )
+  expect_type(bd_mod$description$covariates, "list")
+
+  # check model: the blocks NPex does not define (secondary, initial_conditions, fa) are dropped
+  expect_type(bd_mod$model, "list")
+  expect_equal(
+    names(bd_mod$model),
+    c("primary", "covariates", "lag", "equation", "out", "error")
+  )
+
+  # check support points
+  expect_s3_class(bd_mod$support_point, "data.frame")
+  expect_equal(ncol(bd_mod$support_point), 5)
+  expect_equal(nrow(bd_mod$support_point), 17)
+  expect_equal(names(bd_mod$support_point), c("ka", "ke", "v", "tlag1", "prob"))
+})
+
+# =================================
+# Test checkBDmodel
+# =================================
+
+test_that("checkBDmodel accepts a model built by createBDmodel", {
+  bd_mod <- createBDmodel(
+    PM_result,
+    drug_name = "Amikacin",
+    model_name = "AmikacinPM"
+  )
+
+  expect_true(checkBDmodel(bd_mod))
+})
+
+test_that("checkBDmodel catches what BestDose would refuse", {
+  bd_mod <- createBDmodel(
+    PM_result,
+    drug_name = "Amikacin",
+    model_name = "AmikacinPM"
+  )
+
+  broken <- bd_mod
+  broken$model$primary <- broken$model$primary["ka"]
+  expect_error(checkBDmodel(broken), "At least two primary parameters")
+
+  broken <- bd_mod
+  broken$description$drug <- ""
+  expect_error(checkBDmodel(broken), "A drug must be assigned")
+
+  broken <- bd_mod
+  broken$description$compartment$number <- 3
+  expect_error(
+    checkBDmodel(broken),
+    "does not match the number of named compartments"
+  )
+
+  broken <- bd_mod
+  broken$model$covariates$weight <- NULL
+  expect_error(
+    checkBDmodel(broken),
+    "covariates in the description and in the model block do not match"
+  )
+})
\ No newline at end of file