diff --git a/.Rbuildignore b/.Rbuildignore
index fe62c767..7825d146 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -30,3 +30,5 @@
 ^\.github$
 ^vignettes/articles$
 ^CRAN-SUBMISSION$
+^\.positai$
+^\.claude$
diff --git a/.gitignore b/.gitignore
index ba612b9d..e5de8f87 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,3 +7,4 @@ PopED_output_summary_mfea_opt_1.txt
 inst/doc
 .DS_Store
 docs
+.positai
diff --git a/DESCRIPTION b/DESCRIPTION
index fc320567..d84c4c3f 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -32,7 +32,8 @@ Suggests:
     gridExtra,
     covr,
     devtools,
-    mrgsolve
+    mrgsolve,
+    babelmixr2
 Authors@R: c(
     person("Andrew C.","Hooker", email="andrew.hooker@farmaci.uu.se", 
       role=c("aut","cre","trl","cph"),
@@ -58,7 +59,7 @@ URL: https://andrewhooker.github.io/PopED/, https://github.com/andrewhooker/PopE
 BugReports: https://github.com/andrewhooker/PopED/issues
 Copyright: 2014-2021 Andrew C. Hooker
 Encoding: UTF-8
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 VignetteBuilder: knitr
 Config/Needs/website:
     mrgsolve,
diff --git a/NEWS.md b/NEWS.md
index b783699f..bfc5518b 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,7 @@
 # PopED (development version)
 
+* Make parallelization work with `babelmixr2` models on windows (#79)
+
 # PopED 0.7.0
 
 * `create.poped.database()` now uses a better method of identifying the total number of parameters of each type (bpop, d, sigma, etc.) in a user defined model parameter function (the `ff_fun` argument in `create.poped.database()`) (#73).
diff --git a/R/Doptim.R b/R/Doptim.R
index 85183879..315fd5a2 100644
--- a/R/Doptim.R
+++ b/R/Doptim.R
@@ -74,7 +74,6 @@ Doptim <- function(poped.db,ni, xt, model_switch, x, a, bpopdescr,
                    iter_max=10,
                    ...){
   
-  
   ## update poped.db with options supplied in function
   called_args <- match.call()
   default_args <- formals()
diff --git a/R/LinMatrixH.R b/R/LinMatrixH.R
index 3b241b84..42b36535 100644
--- a/R/LinMatrixH.R
+++ b/R/LinMatrixH.R
@@ -23,7 +23,8 @@
 ## Function translated automatically using 'matlab.to.r()'
 ## Author: Andrew Hooker
 
-LinMatrixH <- function(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,poped.db){
+LinMatrixH <- function(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,
+poped.db){
   #----------Model linearization with respect to epsilon.
   #
   # size of return is (samples per individual x number of epsilons) 
diff --git a/R/LinMatrixL.R b/R/LinMatrixL.R
index 29eba463..cbc19fbc 100644
--- a/R/LinMatrixL.R
+++ b/R/LinMatrixL.R
@@ -22,7 +22,6 @@
 ## Author: Andrew Hooker
 
 LinMatrixL <- function(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,poped.db){
-  
   if((poped.db$parameters$NumRanEff==0)){
     y=0
   } else {
diff --git a/R/LinMatrixLH.R b/R/LinMatrixLH.R
index 27e4e1a5..67f67a28 100644
--- a/R/LinMatrixLH.R
+++ b/R/LinMatrixLH.R
@@ -56,6 +56,7 @@ LinMatrixLH <- function(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,NumEPS,poped
 
 #Helper function to get the hessian for the AD derivative
 new_ferror_file <- function(model_switch,deriv_vec,xt_ind,x,a,bpop,bocc_ind,poped.db){
+
   fg0=feval(poped.db$model$fg_pointer,x,a,bpop,deriv_vec(1:poped.db$parameters$NumRanEff),bocc_ind) #Interaction
   returnArgs <- feval(poped.db$model$ferror_pointer,model_switch,xt_ind,fg0,deriv_vec(poped.db$parameters$NumRanEff+1:length(deriv_vec)),poped.db) 
   f_error <- returnArgs[[1]]
diff --git a/R/RS_opt.R b/R/RS_opt.R
index 703bb182..459b313b 100644
--- a/R/RS_opt.R
+++ b/R/RS_opt.R
@@ -74,7 +74,6 @@ RS_opt <- function(poped.db,
                        compute_inv=TRUE,
                        ...){
   
-  
   # Only get inputs that are needed, not double inputs
   # needed inputs to function: get first then run function
   # poped.db$settings$cfaxt     0.001
diff --git a/R/blockexp.R b/R/blockexp.R
index 4bb6a248..2e6b8b9e 100644
--- a/R/blockexp.R
+++ b/R/blockexp.R
@@ -25,7 +25,7 @@
 blockexp <- function(fn,poped.db,e_flag=FALSE,
                      opt_xt=poped.db$settings$optsw[2],opt_a=poped.db$settings$optsw[4],opt_x=poped.db$settings$optsw[4],
                      opt_samps=poped.db$settings$optsw[1],opt_inds=poped.db$settings$optsw[5]){
-  
+
   fprintf(fn,'==============================================================================\n')
   fprintf(fn,'Model description : %s \n',poped.db$settings$modtit)
   fprintf(fn,'\n')
@@ -155,6 +155,7 @@ blockexp <- function(fn,poped.db,e_flag=FALSE,
 }
 
 print_params <- function (params,name_str, fn, poped.db, param_sqrt=FALSE,head_txt=NULL,matrix_elements=F,e_flag=FALSE) {
+      
   if(is.null(head_txt)) head_txt <- "Parameter Values"
   uncer_txt <- ""
   if(e_flag) uncer_txt <- " (Uncertainty Distribution)"
diff --git a/R/blockopt.R b/R/blockopt.R
index 79338488..386ba513 100644
--- a/R/blockopt.R
+++ b/R/blockopt.R
@@ -20,7 +20,7 @@
 ## Author: Andrew Hooker
 
 blockopt <- function(fn,poped.db,opt_method=""){
-  
+
   if(any(opt_method==c("RS","SG","DO"))){
     fprintf(fn,'==============================================================================\n')
     fprintf(fn,'Optimization Settings\n\n')
@@ -50,4 +50,4 @@ blockopt <- function(fn,poped.db,opt_method=""){
     fprintf(fn,"\n")
   }
   return( ) 
-}
\ No newline at end of file
+}
diff --git a/R/calc_ofv_and_fim.R b/R/calc_ofv_and_fim.R
index 47f2be76..09cba84b 100644
--- a/R/calc_ofv_and_fim.R
+++ b/R/calc_ofv_and_fim.R
@@ -47,7 +47,8 @@ calc_ofv_and_fim <- function (poped.db,
                               ofv_fun = poped.db$settings$ofv_fun,
                               evaluate_fim = TRUE,
                               ...) {
-  
+
+
   ## compute the OFV
   if((ofv==0)){
     if(d_switch){ 
@@ -177,4 +178,4 @@ calc_ofv_and_fim <- function (poped.db,
     fim <- fmf
   }
   return(list(ofv=ofv,fim=fim)) 
-}
\ No newline at end of file
+}
diff --git a/R/create_ofv.R b/R/create_ofv.R
index 03c2cf5b..22163b7c 100644
--- a/R/create_ofv.R
+++ b/R/create_ofv.R
@@ -15,7 +15,6 @@ create_ofv <- function(poped.db,
                        ofv_fun = poped.db$settings$ofv_fun,
                        transform_parameters=T,
                        ...){
-  
   #------------ update poped.db with options supplied in function
   called_args <- match.call()
   default_args <- formals()
diff --git a/R/dfimdalpha.R b/R/dfimdalpha.R
index ee072bfc..b72ad056 100644
--- a/R/dfimdalpha.R
+++ b/R/dfimdalpha.R
@@ -33,4 +33,4 @@ dfimdalpha <- function(alpha, model_switch,groupsize,ni,xtoptn,xoptn,aoptn,bpopd
   }
   return(list( grad= grad,fim =fim )) 
 }
-  
\ No newline at end of file
+  
diff --git a/R/evaluate.e.ofv.fim.R b/R/evaluate.e.ofv.fim.R
index 705ddb39..4ea8d88e 100644
--- a/R/evaluate.e.ofv.fim.R
+++ b/R/evaluate.e.ofv.fim.R
@@ -45,7 +45,6 @@ evaluate.e.ofv.fim <- function(poped.db,
                                use_laplace=poped.db$settings$iEDCalculationType, 
                                laplace.fim=FALSE, 
                                ...){
-
   ## update poped.db with options supplied in function
   called_args <- match.call()
   default_args <- formals()
diff --git a/R/evaluate.fim.R b/R/evaluate.fim.R
index 4b871fde..97c8b7ad 100644
--- a/R/evaluate.fim.R
+++ b/R/evaluate.fim.R
@@ -72,7 +72,6 @@ evaluate.fim <- function(poped.db,
                          deriv.type = NULL,
                          ...){
   
-  
   if(is.null(bpop.val)) bpop.val <- poped.db$parameters$param.pt.val$bpop
   if(is.null(d_full)) d_full <- poped.db$parameters$param.pt.val$d
   if(is.null(docc_full)) docc_full <- poped.db$parameters$param.pt.val$docc
diff --git a/R/evaluate_design.R b/R/evaluate_design.R
index da0826b5..52288d4d 100644
--- a/R/evaluate_design.R
+++ b/R/evaluate_design.R
@@ -23,4 +23,4 @@ evaluate_design <- function(poped.db, ...) {
     colnames(out$fim) <- names(out$rse)
   }
   return(out)
-}
\ No newline at end of file
+}
diff --git a/R/evaluate_fim_map.R b/R/evaluate_fim_map.R
index 4bbdc13f..ac7b33c1 100644
--- a/R/evaluate_fim_map.R
+++ b/R/evaluate_fim_map.R
@@ -37,7 +37,6 @@ evaluate_fim_map <- function(poped.db,
                       num_sim_ids = 1000,
                       use_purrr = FALSE,
                       shrink_mat=F){
-  
   # if (poped.db$design$m > 1) {
   #   warning("Shrinkage should only be computed for a single arm, please adjust your script accordingly.")
   # }
@@ -176,4 +175,4 @@ evaluate_fim_map <- function(poped.db,
   return(out_df)
   
   
-}
\ No newline at end of file
+}
diff --git a/R/evaluate_power.R b/R/evaluate_power.R
index 2925f0d9..06deb4e2 100644
--- a/R/evaluate_power.R
+++ b/R/evaluate_power.R
@@ -85,4 +85,4 @@ evaluate_power <- function(poped.db, bpop_idx, h0=0, alpha=0.05, power=0.80, two
   }
   
   return(out)
-}
\ No newline at end of file
+}
diff --git a/R/get_cv.R b/R/get_cv.R
index 13e0162c..352f5ac0 100644
--- a/R/get_cv.R
+++ b/R/get_cv.R
@@ -62,7 +62,6 @@ get_rse <- function (fim, poped.db,
                      prior_fim = poped.db$settings$prior_fim,
                      #pseudo_on_fail = FALSE,
                      ...) {
-  
   ## update poped.db with options supplied in function
   called_args <- match.call()
   default_args <- formals()
diff --git a/R/get_unfixed_params.R b/R/get_unfixed_params.R
index 813b40df..3a8187f8 100644
--- a/R/get_unfixed_params.R
+++ b/R/get_unfixed_params.R
@@ -15,7 +15,6 @@
 #' @export
 #' @keywords internal
 get_unfixed_params <- function(poped.db,params=NULL){
-  
   if(is.null(params)){
     bpop = poped.db$parameters$bpop[,2,drop=F]
     d = poped.db$parameters$d[,2,drop=F]
diff --git a/R/grad_bpop.R b/R/grad_bpop.R
index 7999b340..cf9d4568 100644
--- a/R/grad_bpop.R
+++ b/R/grad_bpop.R
@@ -13,7 +13,6 @@ grad_bpop <- function(func,select_par,nout,model_switch,xt_ind,x,a,bpop,b_ind,bo
 
 # helper for m2
 helper_v_EBE <- function(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,poped.db) {
-  
   if((poped.db$settings$bCalculateEBE)){
     #zeros(size(b_ind)[1],size(b_ind)[2])
     b_ind_x = ind_estimates(poped.db$mean_data,bpop,d,sigma,t(b_ind),(poped.db$settings$iApproximationMethod==2),FALSE,model_switch,xt_ind,x,a,b_ind,bocc_ind,poped.db)
@@ -28,7 +27,6 @@ helper_v_EBE <- function(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,d,sigma,doc
 
 # helper for m1
 helper_LinMatrix <- function(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,poped.db) {
-  
   epsi0 = zeros(1,length(poped.db$parameters$notfixed_sigma))
   
   # create linearized model
diff --git a/R/graddetmf_ext.R b/R/graddetmf_ext.R
index b8e8190e..1a97209a 100644
--- a/R/graddetmf_ext.R
+++ b/R/graddetmf_ext.R
@@ -2,7 +2,6 @@
 ## Author: Andrew Hooker
 
 graddetmf_ext <- function(model_switch,aX,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db,lndet=FALSE,gradxt=FALSE){
-  
   n = get_fim_size(poped.db)
   m=size(ni,1)
   if (gradxt==FALSE) {
diff --git a/R/gradofv_a.R b/R/gradofv_a.R
index 42d7ae1a..fcdcc402 100644
--- a/R/gradofv_a.R
+++ b/R/gradofv_a.R
@@ -2,7 +2,6 @@
 ## Author: Andrew Hooker
 
 gradofv_a <- function(model_switch,aa,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db){
-  
   #Input: the prior FIM or (empty) and all the other things to calculate the
   #grad with for a
   #Return a vector that is the gradient
diff --git a/R/gradofv_xt.R b/R/gradofv_xt.R
index 2b3f603e..b5ac2687 100644
--- a/R/gradofv_xt.R
+++ b/R/gradofv_xt.R
@@ -2,7 +2,6 @@
 ## Author: Andrew Hooker
 
 gradofv_xt <- function(model_switch,axt,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db){
-  
   #Input: the prior FIM or (empty) and all the other things to calculate the
   #grad with
   #Return a vector that is the gradient and the global structure
diff --git a/R/gradtrmf.R b/R/gradtrmf.R
index 53db2470..77540655 100644
--- a/R/gradtrmf.R
+++ b/R/gradtrmf.R
@@ -3,7 +3,6 @@ gradtrmf <- function(model_switch,aX,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped
   #  Looks at the gradient of tr(FIM^-1) with respect to time (xt) or covariate (a).
   #  problems can arise when a or xt goes negative. So only do forward
   #  differencing. 
-  
   m=size(ni,1)
   if (gradxt == FALSE) {
     gdmf=matrix(1,m,size(a,2))
diff --git a/R/helper.R b/R/helper.R
index e8aa2c01..8623ba47 100644
--- a/R/helper.R
+++ b/R/helper.R
@@ -109,4 +109,4 @@ get_fim_size <- function(poped.db) {
     numnotfixed_covdocc+numnotfixed_sigma+numnotfixed_covsigma
   fim_size <- n_fixed_eff+n_rand_eff
   return(fim_size)
-}
\ No newline at end of file
+}
diff --git a/R/hessian_eta_complex.R b/R/hessian_eta_complex.R
index abd98194..0c813c7d 100644
--- a/R/hessian_eta_complex.R
+++ b/R/hessian_eta_complex.R
@@ -1,6 +1,5 @@
 #Hessian over eta, evaluated at eta_hat => laplace approximation possible
 hessian_eta_complex <- function(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,poped.db,return_gradient=F){
-  
   bAutomatic = FALSE
   epsi0 = zeros(1,length(poped.db$parameters$notfixed_sigma))
   n=length(b_ind)
diff --git a/R/ind_estimates.R b/R/ind_estimates.R
index 3e85e20d..77ac1c73 100644
--- a/R/ind_estimates.R
+++ b/R/ind_estimates.R
@@ -1,7 +1,6 @@
 #Get the emperical bayes estimates for one individual
 #Written for PopED by JN
 ind_estimates <- function(data,bpop,d,sigma,start_bind,bInter,bUDDLike,model_switch,xt_ind,x,a,b_ind,bocc_ind,poped.db){
-  
   b_i = t(start_bind)
   c1 = length(t(start_bind))/2*log(2*pi)
   c2 = 1/2*log(det(d))
diff --git a/R/m3.R b/R/m3.R
index 989da56b..f1b76bc8 100644
--- a/R/m3.R
+++ b/R/m3.R
@@ -4,7 +4,6 @@
 m3 <- function(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,bUseVarSigmaDerivative,poped.db){
   #
   # size: (samps per subject^2 x (number of random effects + number of occasion variances + number of sigmas))
-
   dv_dd = NULL
   dv_covd = NULL
   dv_ddocc = NULL
diff --git a/R/mf3.R b/R/mf3.R
index cc589cca..f4bd1002 100644
--- a/R/mf3.R
+++ b/R/mf3.R
@@ -22,7 +22,6 @@
 ## Author: Andrew Hooker
 
 mf3 <- function(model_switch,xt,x,a,bpop,d,sigma,docc,poped.db){
-
   numnotfixed_bpop = sum(poped.db$parameters$notfixed_bpop)
   numnotfixed_d    = sum(poped.db$parameters$notfixed_d)
   numnotfixed_covd = sum(poped.db$parameters$notfixed_covd)
diff --git a/R/mf7.R b/R/mf7.R
index e7aaa6d5..21cf7acd 100644
--- a/R/mf7.R
+++ b/R/mf7.R
@@ -19,7 +19,6 @@
 ## Author: Andrew Hooker
 
 mf7 <- function(model_switch,xt_ind,x,a,bpop,d,sigma,docc,poped.db){
-  
   #This calculation of FIM divides the calculation up into one calculation
   #per model switch
   
diff --git a/R/mf_all.R b/R/mf_all.R
index 369b6c7d..1793c92d 100644
--- a/R/mf_all.R
+++ b/R/mf_all.R
@@ -3,7 +3,6 @@
 ## Author: Andrew Hooker
 
 mf_all <- function(model_switch,xt_ind,x,a,bpop,d,sigma,docc,poped.db){
-
     returnArgs <- switch(poped.db$settings$iFIMCalculationType+1,
                          mf3(model_switch,xt_ind,x,a,bpop,d,sigma,docc,poped.db), #Default (with no assumption that bpop and b are uncorrelated)
                          mf3(model_switch,xt_ind,x,a,bpop,d,sigma,docc,poped.db), #Reduced FIM
diff --git a/R/mfea.R b/R/mfea.R
index d019f424..0fcb5093 100644
--- a/R/mfea.R
+++ b/R/mfea.R
@@ -36,7 +36,6 @@ mfea <- function(poped.db,model_switch,ni,xt,x,a,bpopdescr,ddescr,maxxt,minxt,ma
                  #opt_inds=poped.db$settings$optsw[5],
                  trflag=T,
                  ...){
-  
   if((poped.db$settings$EACriteria!=1)){
     stop(sprintf('The criteria that can be used is Modified Fedorov Exchange Algorithm (EACriteria = 1)'))
   }
diff --git a/R/ofv_fim.R b/R/ofv_fim.R
index 5e27d4e7..bc52cf6c 100644
--- a/R/ofv_fim.R
+++ b/R/ofv_fim.R
@@ -37,7 +37,6 @@ ofv_fim <- function(fmf,poped.db,
                     ds_index=poped.db$parameters$ds_index,
                     use_log = TRUE,
                     ...){
-  
   #Input: the FIM
   #Return the single value that should be maximized
   
diff --git a/R/optim_ARS.R b/R/optim_ARS.R
index 099930bb..4a059b8a 100644
--- a/R/optim_ARS.R
+++ b/R/optim_ARS.R
@@ -62,6 +62,7 @@ optim_ARS <- function(par,
                       parallel_type=NULL,
                       num_cores = NULL,
                       mrgsolve_model=NULL,
+                      babelmixr2_model=NULL,
                       seed=round(runif(1,0,10000000)),
                       allow_replicates=TRUE,
                       replicates_index=seq(1,length(par)), # same value, parameters can not be the same value
@@ -134,6 +135,7 @@ optim_ARS <- function(par,
                                parallel_type=parallel_type,
                                num_cores=num_cores,
                                mrgsolve_model=mrgsolve_model,
+                               babelmixr2_model=babelmixr2_model,
                                ...) 
     on.exit(if(parallel && (attr(parallel,"type")=="snow")) 
       parallel::stopCluster(attr(parallel,"cluster")))
diff --git a/R/optim_LS.R b/R/optim_LS.R
index 22f7ce69..54ebe371 100644
--- a/R/optim_LS.R
+++ b/R/optim_LS.R
@@ -67,6 +67,7 @@ optim_LS <- function(par,
                      parallel_type=NULL,
                      num_cores = NULL,
                      mrgsolve_model=NULL,
+                     babelmixr2_model=NULL,
                      seed=round(runif(1,0,10000000)),
                      allow_replicates=TRUE,
                      replicates_index=seq(1,length(par)), # same value, parameters can not be the same value
@@ -124,7 +125,7 @@ optim_LS <- function(par,
   
   # start parallel computing
   if(parallel){
-    parallel <- start_parallel(parallel,seed=seed,parallel_type=parallel_type,num_cores=num_cores,mrgsolve_model=mrgsolve_model,...) 
+    parallel <- start_parallel(parallel,seed=seed,parallel_type=parallel_type,num_cores=num_cores,mrgsolve_model=mrgsolve_model,babelmixr2_model=babelmixr2_model,...) 
     on.exit(if(parallel && (attr(parallel,"type")=="snow")) parallel::stopCluster(attr(parallel,"cluster")))
   }   
   #if(is.null(iter_chunk)) if(parallel) iter_chunk <- attr(parallel,"cores") else iter_chunk <- 1
diff --git a/R/optimize_n.R b/R/optimize_n.R
index 403b0bcf..1ae3db57 100644
--- a/R/optimize_n.R
+++ b/R/optimize_n.R
@@ -62,7 +62,6 @@ optimize_groupsize <-
            props = c(poped.db$design$groupsize/sum(poped.db$design$groupsize)),
            trace=1,
            ...){
-    
     # need to fix:
     # limits on the proportions to account for max and min values of N in each group
     # return actual values.
@@ -231,7 +230,6 @@ optimize_n_rse <- function(poped.db,
                            allowed_values = seq(poped.db$design$m,
                                                 sum(poped.db$design$groupsize)*5,
                                                 by=poped.db$design$m)){
-  
   n_per_group = poped.db$design$groupsize
   n_tot <- sum(n_per_group)
   props = c(n_per_group/n_tot)
@@ -275,7 +273,6 @@ optimize_n_eff <- function(poped.db,
                            ofv_ref,
                            norm_group_fim = NULL,  
                            ...){
-  
   n_per_group = poped.db$design$groupsize
   n_tot <- sum(n_per_group)
   props = c(n_per_group/n_tot)
diff --git a/R/par_and_space_tbl.R b/R/par_and_space_tbl.R
index d529ddc6..81f5edb3 100644
--- a/R/par_and_space_tbl.R
+++ b/R/par_and_space_tbl.R
@@ -264,7 +264,7 @@ get_par_and_space_optim <- function(poped.db,
                                     cont_cat = "both",
                                     warn_when_none=TRUE) 
 {
-    
+
   type <- index <- fixed <- cont <- par <- lower <- upper <- NULL
   
   #----------- checks
diff --git a/R/pargen.R b/R/pargen.R
index bde9b221..8374027c 100644
--- a/R/pargen.R
+++ b/R/pargen.R
@@ -29,7 +29,6 @@
 ## Author: Andrew Hooker
 
 pargen <- function(par,user_dist_pointer,sample_size,bLHS,sample_number,poped.db){
-  
   nvar=size(par,1)
   ret=zeros(sample_size,nvar)
   
diff --git a/R/plot_efficiency_of_windows.R b/R/plot_efficiency_of_windows.R
index 64dc21fd..a263929f 100644
--- a/R/plot_efficiency_of_windows.R
+++ b/R/plot_efficiency_of_windows.R
@@ -69,7 +69,6 @@ plot_efficiency_of_windows <- function(poped.db,
                                        #mrgsolve_model=NULL,
                                        seed=round(runif(1,0,10000000)),
                                        ...){
-  
   if(!is.null(seed)) set.seed(seed)
   
   design = poped.db$design
@@ -340,6 +339,7 @@ plot_efficiency_of_windows <- function(poped.db,
                                #parallel_type=parallel_type,
                                #num_cores=num_cores,
                                #mrgsolve_model=mrgsolve_model,
+                               babelmixr2_model=poped.db$babelmixr2,
                                ...) 
     on.exit(if(parallel && (attr(parallel,"type")=="snow")) parallel::stopCluster(attr(parallel,"cluster")))
   }  
@@ -394,4 +394,4 @@ plot_efficiency_of_windows <- function(poped.db,
   p <- p+ylab("Value in %")+xlab("Simulation number of design sampled from defined windows")
   #p <- p + stat_summary(fun.y = "mean", geom="hline")
   return(p)
-}
\ No newline at end of file
+}
diff --git a/R/poped_optim.R b/R/poped_optim.R
index d6e15aff..2ee10d49 100644
--- a/R/poped_optim.R
+++ b/R/poped_optim.R
@@ -94,6 +94,7 @@ poped_optim <- function(poped.db,
                         parallel_type=NULL,
                         num_cores = NULL,
                         mrgsolve_model = NULL,
+                        babelmixr2_model=poped.db$babelmixr2,
                         loop_methods=ifelse(length(method)>1,TRUE,FALSE),
                         iter_max = 10,
                         stop_crit_eff = 1.001,
@@ -106,7 +107,6 @@ poped_optim <- function(poped.db,
                         allow_replicates_a=TRUE,
                         ...){
   
-  
  
   #------------ update argument list with called arguments
   arg.list <- formals()
diff --git a/R/poped_optim_1.R b/R/poped_optim_1.R
index 217bd877..41be6a97 100644
--- a/R/poped_optim_1.R
+++ b/R/poped_optim_1.R
@@ -92,7 +92,6 @@ poped_optim_1 <- function(poped.db,
                         ofv_fun = poped.db$settings$ofv_fun,
                         maximize=T,
                         ...){
-  
   #------------ update poped.db with options supplied in function
   called_args <- match.call()
   default_args <- formals()
diff --git a/R/poped_optim_2.R b/R/poped_optim_2.R
index e77a8422..08895598 100644
--- a/R/poped_optim_2.R
+++ b/R/poped_optim_2.R
@@ -94,7 +94,6 @@ poped_optim_2 <- function(poped.db,
                           maximize=T,
                           transform_parameters = F,
                           ...){
-  
   #------------ update poped.db with options supplied in function
   called_args <- match.call()
   default_args <- formals()
diff --git a/R/poped_optim_3.R b/R/poped_optim_3.R
index 3dee7851..8acf06e1 100644
--- a/R/poped_optim_3.R
+++ b/R/poped_optim_3.R
@@ -96,7 +96,6 @@ poped_optim_3 <- function(poped.db,
                           allow_replicates_xt=TRUE,
                           allow_replicates_a=TRUE,
                           ...){
-  
   #------------ update poped.db with options supplied in function
   called_args <- match.call()
   default_args <- formals()
diff --git a/R/poped_optimize.R b/R/poped_optimize.R
index 7399eef7..117f54e4 100644
--- a/R/poped_optimize.R
+++ b/R/poped_optimize.R
@@ -67,7 +67,6 @@ poped_optimize <- function(poped.db,
                            bLHS=poped.db$settings$bLHS,
                            use_laplace=poped.db$settings$iEDCalculationType,
                            ...){
-  
   ## update poped.db with options supplied in function
   called_args <- match.call()
   default_args <- formals()
diff --git a/R/shrinkage.R b/R/shrinkage.R
index 0eb4e86b..0ee5bbe5 100644
--- a/R/shrinkage.R
+++ b/R/shrinkage.R
@@ -31,7 +31,6 @@ shrinkage <- function(poped.db,
                       use_mc = FALSE,
                       num_sim_ids = 1000,
                       use_purrr = FALSE){
-  
   # if (poped.db$design$m > 1) {
   #   warning("Shrinkage should only be computed for a single arm, please adjust your script accordingly.")
   # }
@@ -200,4 +199,4 @@ shrinkage <- function(poped.db,
   return(out_df)
   
   
-}
\ No newline at end of file
+}
diff --git a/R/start_parallel.R b/R/start_parallel.R
index ca7e5d73..8e665217 100644
--- a/R/start_parallel.R
+++ b/R/start_parallel.R
@@ -32,6 +32,7 @@ start_parallel <- function(parallel=TRUE,
                            seed=NULL,
                            dlls=NULL,
                            mrgsolve_model=NULL,
+                           bablemixr2_model=NULL,
                            #cpp_files=NULL,
                            ...)
 {
@@ -87,6 +88,15 @@ start_parallel <- function(parallel=TRUE,
         }
         parallel::clusterCall(cl, mrgsolve::loadso, x=mrgsolve_model)
       }
+      
+      # load babelmixr2 models in workers using .popedCluster
+      if (!is.null(babelmixr2_model)) {
+        if (!requireNamespace("babelmixr2", quietly=TRUE)) {
+          stop("babelmixr2 package needed for this function to work. Please install it.",
+               call.=FALSE)
+        }
+        parallel::clusterCall(cl, babelmixr2::.popedCluster, babelmixr2_model)
+      }
       # if(!is.null(cpp_files)){
       #   for(i in cpp_files){
       #     parallel::clusterCall(cl, 
@@ -108,4 +118,4 @@ start_parallel <- function(parallel=TRUE,
   }
   
   return(parallel)
-}
\ No newline at end of file
+}
diff --git a/R/v.R b/R/v.R
index 7b0b89b0..f003b972 100644
--- a/R/v.R
+++ b/R/v.R
@@ -3,7 +3,6 @@
 
 v <- function(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,poped.db){
   # number of samples X number of samples (per individual)
-  
   bUseAutoCorrelation = !isempty(poped.db$model$auto_pointer)
   
   bUseFullSigmaCorrelation = FALSE
diff --git a/R/write_iterationfile.R b/R/write_iterationfile.R
index acbc3196..d0759d75 100644
--- a/R/write_iterationfile.R
+++ b/R/write_iterationfile.R
@@ -3,7 +3,6 @@
 ## Author: Andrew Hooker
 
 write_iterationfile <- function(strSearch,iteration,xt,a,x,ni,groupsize,fim,ofv,poped.db){
-
 ## #Make sure that this doesn't crash the search
 ## try 
 ##     fn = file(poped.db$settings$strIterationFileName, 'w')
diff --git a/inst/examples/ex.1.b.PK.1.comp.oral.md.re-parameterize.R b/inst/examples/ex.1.b.PK.1.comp.oral.md.re-parameterize.R
index 22929120..150ea5e0 100644
--- a/inst/examples/ex.1.b.PK.1.comp.oral.md.re-parameterize.R
+++ b/inst/examples/ex.1.b.PK.1.comp.oral.md.re-parameterize.R
@@ -16,12 +16,12 @@ fg.PK.1.comp.oral.md.param.2 <- function(x,a,bpop,b,bocc){
 }
 
 ## -- Define design and design space
-poped.db <- create.poped.database(ff_file="ff.PK.1.comp.oral.md.KE",
-                                  fg_file="fg.PK.1.comp.oral.md.param.2",
-                                  fError_file="feps.add.prop",
+poped.db <- create.poped.database(ff_fun=ff.PK.1.comp.oral.md.KE,
+                                  fg_fun=fg.PK.1.comp.oral.md.param.2,
+                                  fError_fun=feps.add.prop,
                                   groupsize=20,
                                   m=2,
-                                  sigma=c(0.04,5e-6),
+                                  sigma=c(PROP=0.04,ADD=5e-6),
                                   bpop=c(V=72.8,KA=0.25,KE=3.75/72.8,Favail=0.9), 
                                   d=c(V=0.09,KA=0.09,KE=0.25^2), 
                                   notfixed_bpop=c(1,1,1,0),