diff --git a/R/differential_expression.R b/R/differential_expression.R
index e6e9c24..0c9fd5b 100644
--- a/R/differential_expression.R
+++ b/R/differential_expression.R
@@ -93,16 +93,19 @@ compare_expression <- function(x, umi, group, val1, val2, method = 'LRT', bin_si
   for (i in 1:max_bin) {
     genes_bin <- genes[bin_ind == i]
     if (method == 't_test') {
-      bin_res <- future_lapply(genes_bin, function(gene) {
-        model_comparison_ttest(y[gene, use_cells], group)
-      })
+      bin_res <- future_lapply(
+        X = genes_bin, 
+        FUN = function(gene) {model_comparison_ttest(y[gene, use_cells], group)},
+        future.seed = TRUE)
     }
     if (method == 'LRT') {
       mu <- x$model_pars_fit[genes_bin, -1, drop=FALSE] %*% t(regressor_data)  # in log space
       y <- as.matrix(umi[genes_bin, use_cells])
-      bin_res <- future_lapply(genes_bin, function(gene) {
-        model_comparison_lrt(y[gene, ], mu[gene, ], x$model_pars_fit[gene, 'theta'], group, weights)
-      })
+      bin_res <- future_lapply(
+        X = genes_bin, 
+        FUN = function(gene) {
+          model_comparison_lrt(y[gene, ], mu[gene, ], x$model_pars_fit[gene, 'theta'], group, weights)},
+        future.seed = TRUE)
     }
     if (method == 'LRT_reg') {
       LB <- min(x$genes_log_mean_step1)
@@ -172,9 +175,12 @@ compare_expression <- function(x, umi, group, val1, val2, method = 'LRT', bin_si
       y_theta[o] <- 10 ^ ksmooth(x = x$genes_log_mean_step1, y = log10(x$model_pars[, 'theta']),
                                  x.points = y_log_mean, bandwidth = bw, kernel='normal')$y
       names(y_theta) <- genes_bin
-      bin_res <- future_lapply(genes_bin, function(gene) {
-        return(model_comparison_lrt_free3(gene, y[gene, ], y_theta[gene], x$model_str, cell_attr, group, weights, randomize))
-      })
+      bin_res <- future_lapply(
+        X = genes_bin, 
+        FUN = function(gene) {
+          return(model_comparison_lrt_free3(gene, y[gene, ], y_theta[gene], x$model_str, cell_attr, group, weights, randomize))
+        }, 
+        future.seed = TRUE)
     }
     res[[i]] <- do.call(rbind, bin_res)
     if (verbosity > 1) {
diff --git a/R/vst.R b/R/vst.R
index 21198fd..79e7d7f 100644
--- a/R/vst.R
+++ b/R/vst.R
@@ -472,7 +472,8 @@ get_model_pars <- function(genes_step1, bin_size, umi, model_str, cells_step1,
         if (method == "glmGamPoi") {
           return(fit_glmGamPoi(umi = umi_bin_worker, model_str = model_str, data = data_step1))
         }
-      }
+      },
+      future.seed = TRUE
     )
     model_pars[[i]] <- do.call(rbind, par_lst)
 
@@ -500,14 +501,16 @@ get_model_pars_nonreg <- function(genes, bin_size, model_pars_fit, regressor_dat
     genes_bin <- genes[bin_ind == i]
     mu <- tcrossprod(model_pars_fit[genes_bin, -1, drop=FALSE], regressor_data)
     umi_bin <- as.matrix(umi[genes_bin, ])
-    model_pars_nonreg[[i]] <- do.call(rbind,
-                                      future_lapply(genes_bin, function(gene) {
-                                        fam <- negative.binomial(theta = model_pars_fit[gene, 'theta'], link = 'log')
-                                        y <- umi_bin[gene, ]
-                                        offs <- mu[gene, ]
-                                        fit <- glm(as.formula(model_str_nonreg), data = cell_attr, family = fam, offset=offs)
-                                        return(fit$coefficients)
-                                      }))
+    model_pars_nonreg[[i]] <- do.call(
+      rbind,
+      future_lapply(X = genes_bin, 
+                    FUN = function(gene) {
+                      fam <- negative.binomial(theta = model_pars_fit[gene, 'theta'], link = 'log')
+                      y <- umi_bin[gene, ]
+                      offs <- mu[gene, ]
+                      fit <- glm(as.formula(model_str_nonreg), data = cell_attr, family = fam, offset=offs)
+                      return(fit$coefficients)},
+                    future.seed = TRUE))
     if (verbosity > 1) {
       setTxtProgressBar(pb, i)
     }