How to make the results of the maxent model in biomod2 consistent with those in the package dismo in R?

Question

I have prepared a set of point data and a set of raster data (both of which are tutorial data from the Maxent official website). Next, I used the package tidysdm to generate a pseudo absence dataset, and applied this pseudo absence dataset to the Maxent model. At the same time, I set relevant parameters in the model to ensure that the parameters of the maxent model in biomod2 and dismo are consistent. However, in the final results, whether it is the importance of variables, the response curve plot, or the predicted values of raster data for the entire study scope, their results are inconsistent. Here is my code. How should I set up my code to ensure that the maxent models of these two extension packages run exactly the same?

# load packages
if (!require(tidyverse)) install.packages("tidyverse")
if (!require(sf)) install.packages("sf")
if (!require(raster)) install.packages("raster")
if (!require(terra)) install.packages("terra")
if (!require(tidyterra)) install.packages("tidyterra")
if (!require(remotes)) install.packages("remotes")
if (!require(biomod2)) remotes::install_version("biomod2", version = "4.2-3-5")
if (!require(dismo)) install.packages("dismo")
if (!require(tidysdm)) install.packages("tidysdm")

# input presence points
occ <- read.csv("tutorial-data/samples/bradypus.csv") %>%
  dplyr::select(-1)
head(occ)

# input raster file
fnames <- list.files(path = "tutorial-data/layers", pattern = ".asc$", full.names = TRUE)
predictors <- raster::stack(fnames)
predictors

# generate pseudo absence
pts_df <- tidysdm::sample_pseudoabs(data = st_as_sf(occ,coords = c("dd.long","dd.lat"),crs = 4326),
                          raster = terra::subset(rast(predictors),1),
                          n = 3*nrow(occ),
                          method = "random")

# get presence points
occ <- pts_df %>% 
  dplyr::filter(class == "presence") %>% 
  st_coordinates() %>% 
  as.data.frame() %>% 
  setNames(c("dd.long","dd.lat"))
# get pseudo absence points
abs <- pts_df %>% 
  dplyr::filter(class == "pseudoabs") %>% 
  st_coordinates() %>% 
  as.data.frame() %>% 
  setNames(c("dd.long","dd.lat"))

# generate maxent model with package "dismo"
# Here set the explain variable "ecoreg" into factor type
maxent_arg <- c("betamultiplier=1", "noautofeature", "noquadratic", "noproduct", "nothreshold", "nohinge")
me <- dismo::maxent(
  x = predictors,
  p = occ,
  a = abs,
  factors = "ecoreg",
  removeDuplicates = FALSE,
  path = "testrun",
  args = maxent_arg
)
# make prediction on raster file
r <- predict(me, predictors)


# generate biomod2 presence and absence points data
biomod2_df <- occ %>% 
  mutate(PA1 = TRUE) %>% 
  mutate(resp.var = 1) %>% 
  rbind(abs %>% 
          mutate(PA1 = TRUE) %>% 
          mutate(resp.var = NA)) %>% 
  rename(lon = dd.long,lat = dd.lat) %>% 
  as.data.frame()

# input raster file and transform the variable ecoreg into factor
predictors <- list.files(path = "tutorial-data/layers", pattern = ".asc$", full.names = TRUE) %>% 
  rast() %>% 
  mutate(across(.cols = ecoreg,.fns = as.factor))

myBiomodData <- BIOMOD_FormatingData(resp.var = biomod2_df$resp.var,
                     expl.var = predictors,
                     resp.xy = biomod2_df[,c("lon","lat")],
                     resp.name = "occ",
                     PA.strategy = 'user.defined',
                     PA.user.table = dplyr::select(biomod2_df,PA1))

mod <- "MAXENT"
maxentlqpht <- c(T, F, F, F, F)
betamult <- 1
myBiomodOption <- BIOMOD_ModelingOptions(
  MAXENT = list(
    path_to_maxent.jar = "maxent.jar",
    memory_allocated = NULL,
    maximumiterations = 500,
    betamultiplier = betamult,
    beta_lqp = 0.05,
    linear = maxentlqpht[1],
    quadratic = maxentlqpht[2],
    product = maxentlqpht[3],
    hinge = maxentlqpht[4],
    threshold = maxentlqpht[5]
  )
)
myBiomodModelOut <- BIOMOD_Modeling(
  bm.format = myBiomodData,
  models = mod,
  bm.options = myBiomodOption,
  var.import = 5,
  metric.eval = c("KAPPA", "TSS", "ROC"),
  do.full.models = F,
  CV.perc = .25,
  modeling.id = "biomod2_test"
)

myMODRespPlot2D <-
  bm_PlotResponseCurves(
    bm.out = myBiomodModelOut,
    models.chosen = get_built_models(myBiomodModelOut, run = "RUN1"),
    fixed.var = "median",
    Data = get_formal_data(myBiomodModelOut, "expl.var"),
    show.variables = get_formal_data(myBiomodModelOut, "expl.var.names"),
    do.bivariate = FALSE
  )

response(me)

plot(me)

get_variables_importance(myBiomodModelOut) %>% 
  group_by(expl.var) %>% 
  summarise(var.imp = mean(var.imp)) %>% 
  arrange(desc(var.imp))

myBiomodProj <- BIOMOD_Projection(
  bm.mod = myBiomodModelOut,
  proj.name = "maxent_biomod2",
  new.env = predictors,
  models.chosen = "all",
  metric.binary = "all",
  metric.filter = "all",
  build.clamping.mask = TRUE,
  output.format = ".img",
  do.stack = FALSE
)


temp_r <- rast(c(
  biomod2 = rast("proj_maxent_biomod2_occ_PA1_RUN1_MAXENT.img"),
  dismo = app(x = rast(r), fun = function(x) ifelse(is.na(x), x, 1000 * x))
)) %>%
  mutate(difference = biomod2 - dismo) %>%
  mutate(across(.cols = everything(), .fns = as.numeric))
temp_r
ggplot() +
  geom_spatraster(data = temp_r) +
  facet_wrap(~lyr) +
  scale_fill_whitebox_c(
    palette = "muted"
  ) +
  theme_light()