diff --git a/R/sufficient_stats_methods.R b/R/sufficient_stats_methods.R
index f200404c..c5078193 100644
--- a/R/sufficient_stats_methods.R
+++ b/R/sufficient_stats_methods.R
@@ -332,22 +332,41 @@ update_variance_components.ss <- function(data, params, model, ...) {
     omega     <- matrix(rep(omega_res$diagXtOmegaX, L), nrow = L, ncol = data$p, byrow = TRUE) +
       matrix(rep(1 / model$V, data$p), nrow = L, ncol = data$p, byrow = FALSE)
 
-    # Update the sparse effect variance
-    sparse_var <- mean(colSums(model$alpha * model$V))
-    
+    # Update sigma2 and tau2 via MoM
+    mom_result <- mom_unmappable(data, params, model, omega, model$tau2)
+
+    # If tau2 = 0 skip ash entirely
+    if (mom_result$tau2 == 0) {
+      return(list(
+        sigma2 = mom_result$sigma2,
+        tau2   = 0,
+        theta  = rep(0, data$p)
+      ))
+    }
+
     # Remove the sparse effects
     b <- colSums(model$alpha * model$mu)
     residuals <- data$y - data$X %*% b
 
-    # Build ASH grid based on sparse variance and scaling factor
-    if (sparse_var * params$ash_scaling_factor > 1e-8) {
-      grid_factors <- exp(seq(log(0.01), log(10), length.out = 20 - 1))
-      est_sa2 <- c(0, sparse_var * params$ash_scaling_factor * grid_factors)
+    # Get minimum of active sparse effect variances as upper bound
+    active_V <- model$V[model$V > 0.01]
+    if (length(active_V) > 0) {
+      min_sparse_var <- min(active_V)
     } else {
-      est_sa2 <- c(0, (2^(0.05*(1:20-1)) - 1)^4)
-      est_sa2 <- est_sa2 * (0.01 / max(est_sa2))
+      min_sparse_var <- 0.01  # Default when no active sparse effects
     }
 
+    # Region 1: Very dense near zero (tau2/100 to tau2/10)
+    small_grid <- exp(seq(log(mom_result$tau2/100), log(mom_result$tau2/10), length.out = 10))
+
+    # Region 2: Moderate density around tau2 (tau2/10 to tau2*3)
+    medium_grid <- exp(seq(log(mom_result$tau2/10), log(mom_result$tau2*3), length.out = 6))
+
+    # Region 3: Sparse coverage towards minimum sparse effect
+    large_grid <- exp(seq(log(mom_result$tau2*3), log(min_sparse_var), length.out = 3))
+
+    est_sa2 <- c(0, unique(c(small_grid, medium_grid, large_grid)))
+
     # Call mr.ash with residuals
     mrash_output <- mr.ash.alpha::mr.ash(
       X             = data$X,
@@ -355,7 +374,7 @@ update_variance_components.ss <- function(data, params, model, ...) {
       sa2           = est_sa2,
       intercept     = FALSE,
       standardize   = FALSE,
-      sigma2        = model$sigma2,
+      sigma2        = mom_result$sigma2,
       update.sigma2 = params$estimate_residual_variance,
       max.iter      = 3000
     )
diff --git a/R/susie.R b/R/susie.R
index 1abef5e5..e30aecf9 100644
--- a/R/susie.R
+++ b/R/susie.R
@@ -279,7 +279,6 @@ susie <- function(X, y, L = min(10, ncol(X)),
                   estimate_prior_variance = TRUE,
                   estimate_prior_method = c("optim", "EM", "simple"),
                   unmappable_effects = c("none", "inf", "ash"),
-                  ash_scaling_factor = 0.15,
                   check_null_threshold = 0,
                   prior_tol = 1e-9,
                   residual_variance_upperbound = Inf,
@@ -311,7 +310,7 @@ susie <- function(X, y, L = min(10, ncol(X)),
     prior_weights, null_weight, standardize, intercept,
     estimate_residual_variance, estimate_residual_method,
     estimate_prior_variance, estimate_prior_method,
-    unmappable_effects, ash_scaling_factor, check_null_threshold, prior_tol,
+    unmappable_effects, check_null_threshold, prior_tol,
     residual_variance_upperbound, model_init, coverage,
     min_abs_corr, compute_univariate_zscore, na.rm,
     max_iter, tol, convergence_method, verbose, track_fit,
diff --git a/R/susie_constructors.R b/R/susie_constructors.R
index 78ae0f22..355caa13 100644
--- a/R/susie_constructors.R
+++ b/R/susie_constructors.R
@@ -24,7 +24,6 @@ individual_data_constructor <- function(X, y, L = min(10, ncol(X)),
                                         estimate_prior_variance = TRUE,
                                         estimate_prior_method = "optim",
                                         unmappable_effects = "none",
-                                        ash_scaling_factor = 0.15,
                                         check_null_threshold = 0,
                                         prior_tol = 1e-9,
                                         residual_variance_upperbound = Inf,
@@ -141,7 +140,6 @@ individual_data_constructor <- function(X, y, L = min(10, ncol(X)),
     estimate_prior_variance = estimate_prior_variance,
     estimate_prior_method = estimate_prior_method,
     unmappable_effects = unmappable_effects,
-    ash_scaling_factor = ash_scaling_factor,
     check_null_threshold = check_null_threshold,
     prior_tol = prior_tol,
     residual_variance_upperbound = residual_variance_upperbound,
diff --git a/README.md b/README.md
index 00a59768..6bf6e017 100644
--- a/README.md
+++ b/README.md
@@ -85,11 +85,11 @@ If you use any of the summary data methods such as `susie_ss` or
 If you use the Servin-Stephens prior on residual variance estimates
 (`estimate_residual_method = "Servin_Stephens"`), please also cite:
 
-> Denault, W.R.P., Carbonetto, P., Li, R., Consortium, A.D.F.G., Wang,
-> G. & Stephens, M. (2025). Accounting for uncertainty in residual
-> variances improves calibration of the "Sum of Single Effects" model
-> for small sample sizes. *bioRxiv*, 2025-05. Under review for *Nature
-> Methods*.
+> Denault, W.R.P., Carbonetto, P., Li, R., Alzheimer's Disease Functional
+> Genomics Consortium, Wang, G. & Stephens, M. (2025). Accounting for
+> uncertainty in residual variances improves calibration of the "Sum of
+> Single Effects" model for small sample sizes. *bioRxiv*, 2025-05.
+> Under review for *Nature Methods*.
 
 If you use infinitesimal effects modeling (`unmappable_effects = "inf"`), 
 please also cite:
diff --git a/inst/CITATION b/inst/CITATION
index c552176b..4e1d1e16 100644
--- a/inst/CITATION
+++ b/inst/CITATION
@@ -68,18 +68,18 @@ bibentry(header = paste("If estimate_residual_method = 'Servin_Stephens' is used
          author = c(person("William", "R.P. Denault"),
                     person("Peter", "Carbonetto"),
                     person("Ruixi", "Li"),
-                    person("A.D.F.G.", "Consortium"),
+                    person("Alzheimer's Disease Functional Genomics Consortium"),
                     person("Gao", "Wang"),
                     person("Matthew", "Stephens")),
          journal = "bioRxiv",
          year    = "2025",
          note    = "Under review for Nature Methods",
          textVersion =
-           paste("Denault, W.R.P., Carbonetto, P., Li, R., Consortium, A.D.F.G.,",
-                 "Wang, G. & Stephens, M. (2025). Accounting for uncertainty in",
-                 "residual variances improves calibration of the 'Sum of Single",
-                 "Effects' model for small sample sizes. bioRxiv, 2025-05.",
-                 "Under review for Nature Methods."))
+           paste("Denault, W.R.P., Carbonetto, P., Li, R., Alzheimer's Disease",
+                 "Functional Genomics Consortium, Wang, G. & Stephens, M. (2025).",
+                 "Accounting for uncertainty in residual variances improves",
+                 "calibration of the 'Sum of Single Effects' model for small",
+                 "sample sizes. bioRxiv, 2025-05. Under review for Nature Methods."))
 
 bibentry(header = paste("If unmappable_effects = 'inf' is used,",
                         "please also cite:"),