diff --git a/nimbleModel/R/MCMC_configuration.R b/nimbleModel/R/MCMC_configuration.R
index 9612b9a..ff887fc 100644
--- a/nimbleModel/R/MCMC_configuration.R
+++ b/nimbleModel/R/MCMC_configuration.R
@@ -28,7 +28,7 @@ samplerConfClass <- R6Class(
       if (name == "crossLevel") {
         control <<- c(
           control,
-          list(dependent_nodes = getNodes(model, getDependencies(model$modelDef, target, self = FALSE, nodesAsChars = FALSE), stochOnly = TRUE, nodesAsChars = FALSE))
+          list(dependent_nodes = getNodes(model, getDependencies(model, target, self = FALSE, nodesAsChars = FALSE), stochOnly = TRUE, nodesAsChars = FALSE))
         )
       } # special case for printing dependents of crossLevel sampler (only)
     },
diff --git a/nimbleModel/R/MCMC_conjugacy.R b/nimbleModel/R/MCMC_conjugacy.R
index 67da35d..3f2ecb2 100644
--- a/nimbleModel/R/MCMC_conjugacy.R
+++ b/nimbleModel/R/MCMC_conjugacy.R
@@ -239,7 +239,7 @@ conjugacyRelationshipsClass <- R6Class(
 
       # CHECK: is this sufficiently efficient? Is getting a char representation the best strategy?
       targetNode <- getNodes(model, nodes = nodeRange$toNodeChars(1), nodesAsChars = FALSE)[[1]] # First node as representative.
-      deps <- getNodes(model, getDependencies(model$modelDef, targetNode, self = FALSE, nodesAsChars = FALSE), stochOnly = TRUE, nodesAsChars = FALSE)
+      deps <- getNodes(model, getDependencies(model, targetNode, self = FALSE, nodesAsChars = FALSE), stochOnly = TRUE, nodesAsChars = FALSE)
 
       depTypes <- sapply(deps, function(x) conjugacyObj$checkConjugacyOneDep(model, targetNode, x, restrictLink))
 
@@ -396,8 +396,8 @@ conjugacyClass <- R6Class(
     genSetupFunction = function(dependentCounts, doDependentScreen = FALSE) {
       functionBody <- codeBlockClass()
       functionBody$addCode({
-        calcNodes <- getNodes(model, getDependencies(model$modelDef, target, nodesAsChars = FALSE), nodesAsChars = FALSE)
-        calcNodesDeterm <- getNodes(model, getDependencies(model$modelDef, target, nodesAsChars = FALSE), determOnly = TRUE, nodesAsChars = FALSE)
+        calcNodes <- getNodes(model, getDependencies(model, target, nodesAsChars = FALSE), nodesAsChars = FALSE)
+        calcNodesDeterm <- getNodes(model, getDependencies(model, target, nodesAsChars = FALSE), determOnly = TRUE, nodesAsChars = FALSE)
       })
 
       # if this conjugate sampler is for a multivariate node (i.e., nDim > 0), then we need to determine the size (d)
diff --git a/nimbleModel/R/dataRules.R b/nimbleModel/R/dataRules.R
index 0e1e23e..e294130 100644
--- a/nimbleModel/R/dataRules.R
+++ b/nimbleModel/R/dataRules.R
@@ -151,7 +151,7 @@ makeRulePieces <- function(elements, varName, all, sequenceThreshold = 0.1) {
 }
 
 
-excludeFromPredictiveRules <- function(modelDef, currentRanges, candidateRules) {
+excludeFromPredictiveRules <- function(model, currentRanges, candidateRules) {
   if (!length(candidateRules)) {
     return(NULL)
   }
@@ -164,8 +164,8 @@ excludeFromPredictiveRules <- function(modelDef, currentRanges, candidateRules)
     } else {
       candidateRules[[varName]] <- NULL
     }
-    parents <- getParents(modelDef, range, nodesAsChars = FALSE)
-    candidateRules <- excludeFromPredictiveRules(modelDef, parents, candidateRules)
+    parents <- getParents(model, range, nodesAsChars = FALSE)
+    candidateRules <- excludeFromPredictiveRules(model, parents, candidateRules)
   }
   return(candidateRules)
 }
diff --git a/nimbleModel/R/graphRules.R b/nimbleModel/R/graphRules.R
index 4140204..2bfa84b 100644
--- a/nimbleModel/R/graphRules.R
+++ b/nimbleModel/R/graphRules.R
@@ -98,6 +98,7 @@ graphRuleClass <- R6Class(
       applyGraphRule(fromVarRange, self, removeDuplicates = removeDuplicates)
     },
     getFromRange = function() {
+      # Returns maximal extent; so for y[2:5], it returns 1:5; same for y[c(2,5)].
       if (!length(indexSets$fromIndexSlotToSet)) { # no indexing
         varRange <- varRangeClass$new(fromVarName)
       } else {
diff --git a/nimbleModel/R/modelBaseClass.R b/nimbleModel/R/modelBaseClass.R
index 23f794e..ee85b73 100644
--- a/nimbleModel/R/modelBaseClass.R
+++ b/nimbleModel/R/modelBaseClass.R
@@ -2,6 +2,7 @@
 modelBase_nClass <- nClass(
   classname = "modelBase_nClass",
   Rpublic = list(
+    thisModel = NULL,
     modelDef = NULL,
     dataRules = NULL,
     nondataRules = NULL,
@@ -13,6 +14,8 @@ modelBase_nClass <- nClass(
       if (isTRUE(.GlobalEnv$.debugModelInit)) browser()
       super$initialize(...)
 
+      self$thisModel <- self   # A hack needed because we use `self` as method argument name below.
+
       # TODO: is there a better way to populate declFunNameToIndex in Cpublic?
       declFunNameToIndex <- self$declFunNameToIndex_
 
@@ -111,7 +114,7 @@ modelBase_nClass <- nClass(
           dataRule$rule$apply(dataRule$varName)
         })
       }))
-      self$predictiveRules <- excludeFromPredictiveRules(modelDef, dataRanges, candidateRules)
+      self$predictiveRules <- excludeFromPredictiveRules(self, dataRanges, candidateRules)
 
       # nonpredictive rules
       candidateRules <- unlist(lapply(modelDef$calcRules, function(oneVarRules) {
@@ -365,24 +368,36 @@ modelBase_nClass <- nClass(
         return(expr)
       }
     },
-    getDependencies = function(nodes, self = TRUE, downstream = FALSE, immediateOnly = FALSE,
+    # `self` arg masks the reference to the object.
+    # TODO: perhaps we should rename the arg `includeSelf`, but that is not back compatible.
+    getDependencies = function(nodes, self = TRUE, determOnly = FALSE, stochOnly = FALSE,
+                               includeData = TRUE, dataOnly = FALSE,
+                               includePredictive = nimble::getNimbleOption('getDependenciesIncludesPredictiveNodes'),
+                               predictiveOnly = FALSE, includeRHSonly = FALSE,
+                               downstream = FALSE, immediateOnly = FALSE,
                                nodesAsChars = getNimbleModelOption("nodesAsChars"),
                                returnScalarComponents = FALSE, .sort = FALSE) {
       nimbleModel::getDependencies(
-        modelDef, nodes, self, downstream, immediateOnly,
+        thisModel, nodes, self,
+        determOnly, stochOnly, includeData, dataOnly,
+        includePredictive, predictiveOnly, includeRHSonly,
+        downstream, immediateOnly,
         nodesAsChars, returnScalarComponents, .sort
       )
     },
-    getParents = function(nodes, self = FALSE, upstream = FALSE, immediateOnly = FALSE,
+    getParents = function(nodes, self = FALSE,
+                          determOnly = FALSE, stochOnly = FALSE,
+                          includeData = TRUE, dataOnly = FALSE, includeRHSonly = FALSE,
+                          upstream = FALSE, immediateOnly = FALSE,
                           nodesAsChars = getNimbleModelOption("nodesAsChars"),
                           returnScalarComponents = FALSE, .sort = FALSE) {
       nimbleModel::getParents(
-        modelDef, nodes, self, upstream, immediateOnly,
+        thisModel, nodes, self,
+        determOnly, stochOnly, includeData, dataOnly, includeRHSonly,        
+        upstream, immediateOnly,
         nodesAsChars, returnScalarComponents, .sort
       )
     },
-    # TODO: not working because `nimbleModel::getNodes` needs the model not just modelDef.
-    # Once we integrate modelClass with modelBase_nClass, we should be able to pass `self`.
     getNodes = function(nodes, determOnly = FALSE, stochOnly = FALSE,
                         includeData = TRUE, dataOnly = FALSE,
                         includeRHSonly = FALSE,
diff --git a/nimbleModel/R/modelDef.R b/nimbleModel/R/modelDef.R
index 30b6390..590893c 100644
--- a/nimbleModel/R/modelDef.R
+++ b/nimbleModel/R/modelDef.R
@@ -1024,52 +1024,6 @@ modelDefClass <- R6Class(
 )
 
 
-# Core graph and node querying functions in the model API.
-# These are standalone functions for now, but may become
-# part of model class. That said, more naturally part of modelDef class.
-
-# TODO: move these functions into a new stand-alone code file for user-facing functions?
-
-# Note: `getDependencies` and `getParents` cannot handle `stochOnly` or `determOnly`
-# because a given varRange result for getParents could be partially stochastic and
-# partially deterministic. Instead a user would pass the result through `getNodes()`.
-# Similarly, filtering by RHSonly will be done in `getNodes()`.
-
-# Note: data-related flags not handled as that relates to flags on a model
-# and not part of modelDef.
-
-# TODO: these should presumably take the model not modelDef as the first arg.
-# Once we integrate modelClass with modelBase_nClass, we should be able to
-# pass `self` from the getDeps and getParents methods to these functions.
-
-getDependencies <- function(modelDef, nodes,
-                            self = TRUE,
-                            downstream = FALSE, immediateOnly = FALSE,
-                            nodesAsChars = getNimbleModelOption("nodesAsChars"),
-                            returnScalarComponents = FALSE, .sort = FALSE) {
-  traverseGraph(modelDef$downstreamRules, modelDef$declRules,
-    nodes = nodes,
-    down = TRUE, self = self,
-    follow = downstream, immediateOnly = immediateOnly,
-    nodesAsChars = nodesAsChars, returnScalarComponents = returnScalarComponents,
-    .sort = .sort, modelDef = modelDef
-  )
-}
-
-getParents <- function(modelDef, nodes,
-                       self = FALSE,
-                       upstream = FALSE, immediateOnly = FALSE,
-                       nodesAsChars = getNimbleModelOption("nodesAsChars"),
-                       returnScalarComponents = FALSE, .sort = FALSE) {
-  traverseGraph(modelDef$upstreamRules, modelDef$declRules,
-    nodes = nodes,
-    down = FALSE, self = self,
-    follow = upstream, immediateOnly = immediateOnly,
-    nodesAsChars = nodesAsChars, returnScalarComponents = returnScalarComponents,
-    .sort = .sort, modelDef = modelDef
-  )
-}
-
 
 # Evaluates `if` statements in model code to generate actual model code
 # without any `if` statements. Condition of if statement can use variables
diff --git a/nimbleModel/R/modelFunctions.R b/nimbleModel/R/modelFunctions.R
index ef59532..57b87cf 100644
--- a/nimbleModel/R/modelFunctions.R
+++ b/nimbleModel/R/modelFunctions.R
@@ -216,6 +216,46 @@ taggedClass <- R6Class(
   )
 )
 
+getDependencies <- function(model, nodes,
+                            self = TRUE, determOnly = FALSE, stochOnly = FALSE,
+                            includeData = TRUE, dataOnly = FALSE,
+                            includePredictive = nimble::getNimbleOption('getDependenciesIncludesPredictiveNodes'),
+                            predictiveOnly = FALSE, includeRHSonly = FALSE,
+                            downstream = FALSE, immediateOnly = FALSE,
+                            nodesAsChars = getNimbleModelOption("nodesAsChars"),
+                            returnScalarComponents = FALSE, .sort = FALSE) {
+  traverseGraph(model$modelDef$downstreamRules, model$modelDef$declRules,
+    nodes = nodes,
+    down = TRUE, self = self,
+    determOnly = determOnly, stochOnly = stochOnly,
+    includeData = includeData, dataOnly = dataOnly, includePredictive = includePredictive,
+    predictiveOnly = predictiveOnly, includeRHSonly = includeRHSonly,
+    follow = downstream, immediateOnly = immediateOnly,
+    nodesAsChars = nodesAsChars, returnScalarComponents = returnScalarComponents,
+    .sort = .sort, model = model
+  )
+}
+
+getParents <- function(model, nodes,
+                       self = FALSE, determOnly = FALSE, stochOnly = FALSE,
+                       includeData = TRUE, dataOnly = FALSE, includeRHSonly = FALSE,
+                       upstream = FALSE, immediateOnly = FALSE,
+                       nodesAsChars = getNimbleModelOption("nodesAsChars"),
+                       returnScalarComponents = FALSE, .sort = FALSE) {
+  traverseGraph(model$modelDef$upstreamRules, model$modelDef$declRules,
+    nodes = nodes,
+    down = FALSE, self = self,
+    determOnly = determOnly, stochOnly = stochOnly,
+    includeData = includeData, dataOnly = dataOnly, includePredictive = TRUE,
+    predictiveOnly = FALSE, includeRHSonly = includeRHSonly,
+    follow = upstream, immediateOnly = immediateOnly,
+    nodesAsChars = nodesAsChars, returnScalarComponents = returnScalarComponents,
+    .sort = .sort, model = model
+  )
+}
+
+
+
 # This may not optimally aggregate in cases without contiguity - e.g., 2:4 + 6:8 will become a matrix, even though
 # it may be more efficient to leave it as two nodeRanges.
 #' @export
@@ -236,6 +276,7 @@ aggregate_nodes <- function(nodeSet) {
   names(nodeSet) <- declIDs
   nodeIDs <- lapply(nodeSet, \(x) x$getIDs())
   IDsByDecl <- lapply(split(nodeIDs, declIDs), \(x) unique(nimbleModel:::flatten(x)))
+  IDsByDecl <- IDsByDecl[unique(declIDs)]  # Try to keep in order provided.
   nms <- names(IDsByDecl)
   newNodeSet <- lapply(seq_along(IDsByDecl), \(i) {
     if(sum(nms[i] == declIDs) > 1) {
@@ -288,6 +329,16 @@ setdiff_nodes <- function(nodeSet1, nodeSet2) {
   }
   declIDs1 <- sapply(nodeSet1, \(x) x$decl$declRule$ID)
   declIDs2 <- sapply(nodeSet2, \(x) x$decl$declRule$ID)
+
+  nullCases <- sapply(declIDs1, is.null)
+  RHSonly <- nodeSet1[nullCases]
+  nodeSet1 <- nodeSet1[!nullCases]
+  declIDs1 <- unlist(declIDs1[!nullCases])
+
+  nullCases <- sapply(declIDs2, is.null)
+  nodeSet2 <- nodeSet2[!nullCases]
+  declIDs2 <- unlist(declIDs2[!nullCases])
+
   nodeIDs1 <- lapply(nodeSet1, \(x) x$getIDs())
   nodeIDs2 <- lapply(nodeSet2, \(x) x$getIDs())
   excludeNodeIDs <- lapply(split(nodeIDs2, declIDs2), \(x) unique(nimbleModel:::flatten(x)))
@@ -301,7 +352,7 @@ setdiff_nodes <- function(nodeSet1, nodeSet2) {
       )
     }
   }
-  return(newNodeSet1[!sapply(newNodeSet1, is.null)])
+  return(c(newNodeSet1[!sapply(newNodeSet1, is.null)], RHSonly))
 }
 
 #' @export
@@ -922,8 +973,7 @@ splitLatents <- function(model, paramNodes, latentNodes, calcNodes, calcNodesOth
     latentNodes <- margNodes$randomEffectsNodes
     deps <- model$getNodes(model$getDependencies(latentNodes, self = FALSE, nodesAsChars = FALSE), includeData = FALSE, nodesAsChars = FALSE)
     ## By default, we treat "siblings" of latent nodes as latents.
-    ## This attempts to have fixed effects in latents,
-    ## along with random effects.
+    ## This attempts to have fixed effects in latents, along with random effects.
     newLatents <- model$getNodes(model$getParents(deps, nodesAsChars = FALSE), stochOnly = TRUE, includeData = FALSE, nodesAsChars = FALSE)
     paramNodes <- setdiff_nodes(paramNodes, newLatents)
     latentNodes <- aggregate_nodes(c(latentNodes, newLatents))
diff --git a/nimbleModel/R/processModelGraph.R b/nimbleModel/R/processModelGraph.R
index 076ecc2..49a01b6 100644
--- a/nimbleModel/R/processModelGraph.R
+++ b/nimbleModel/R/processModelGraph.R
@@ -236,14 +236,15 @@ setSortIDs <- function(calcRules) {
 # By default stops at stochastic nodes, unless requested to go through
 # (`follow = TRUE`) or to stop at immediate parent or child
 # (`immediateOnly = TRUE`).
-# Result is a set of varRanges (not nodeRanges), so users may need to
-# pass result through `getNodes`.
 traverseGraph <- function(streamRules, declRules,
                           nodes, down, self = TRUE,
+                          determOnly = FALSE, stochOnly = FALSE,
+                          includeData = TRUE, dataOnly = FALSE,
+                          includePredictive = TRUE, predictiveOnly = FALSE, includeRHSonly = FALSE,
                           follow = FALSE, immediateOnly = FALSE,
                           nodesAsChars = getNimbleModelOption("nodesAsChars"),
                           returnScalarComponents = FALSE, .sort = FALSE,
-                          modelDef = NULL) {
+                          model = NULL) {
   # A single varRange can have elements that don't share a sortID when converted to calcRange representation,
   # so we can't sort varRanges.
   if (.sort && !nodesAsChars) {
@@ -254,129 +255,57 @@ traverseGraph <- function(streamRules, declRules,
 
   results <- traverseGraphRecurse(streamRules, nodes, down, follow, immediateOnly)
 
-  # Need to handle "self" for three cases: (a) when an input node is a full variable,
-  # (b) character expression for a range, or (c) an actual varRange or nodeRange.
-
-  varNames <- unlist(lapply(nodes, getVarName))
-  vars <- nodes == varNames
-  selfRangeFromVars <- flatten(lapply(
-    nodes[vars],
-    function(varName) {
-      lapply(
-        declRules[[varName]]$rules,
-        function(declRule) declRule$fullRange
-      )
-    }
-  ))
-
-  # CHECK: can we just create the varRange from the char string and not use `declRule$apply`
-  # to make sure we have decl-specific varRanges?
-  charRanges <- is.character(nodes) & !vars
-  selfRangeFromCharRanges <- flatten(lapply(
-    nodes[charRanges],
-    function(node) {
-      lapply(
-        declRules[[getVarName(node)]]$rules,
-        function(declRule) {
-          tmp <- declRule$apply(node)
-          if (is.null(tmp)) NULL else tmp$toVarRange()
-        }
-      )
-    }
-  ))
-  if (identical(selfRangeFromCharRanges, list(NULL))) {
-    selfRangeFromCharRanges <- NULL
-  }
-  selfRangeFromNodes <- lapply(
-    nodes[!vars & !charRanges],
-    function(node) {
-      if (inherits(node, "nodeRangeClass")) {
-        return(node$toVarRange())
-      } else {
-        return(node)
-      }
-    }
-  )
-  selfRanges <- c(selfRangeFromNodes, selfRangeFromVars, selfRangeFromCharRanges)
-
-  if (length(results)) {
-    # Exclude self (add back below if needed).
-    # This helps avoid duplication, though that might be handled fully by removeDuplicateVarRanges.
-    for (i in seq_along(selfRanges)) {
-      resultsNames <- unlist(lapply(results, function(x) x$varName))
-      wh <- which(resultsNames == selfRanges[[i]]$varName)
-      if (length(wh)) {
-        newResults <- list()
-        for (idx in wh) {
-          newResults <- c(
-            newResults,
-            lapply(
-              exclude(results[[idx]], selfRanges[[i]]),
-              function(rule) rule$fullRange
-            )
-          )
-        }
-        results <- c(results[-wh], newResults)
-      }
-    }
-  }
+  results <- model$getNodes(results, determOnly = determOnly, stochOnly = stochOnly,
+                            includeData = includeData, dataOnly = dataOnly,
+                            includePredictive = includePredictive, predictiveOnly = predictiveOnly,
+                            includeRHSonly = includeRHSonly, nodesAsChars = FALSE)
 
   if (self) {
-    results <- c(selfRanges, results)
-  }
-
-  if (!length(results)) {
-    return(NULL)
+    results <- aggregate_nodes(c(model$getNodes(nodes, nodesAsChars = FALSE), results))
+  } else {
+    results <- setdiff_nodes(aggregate_nodes(results), model$getNodes(nodes, nodesAsChars = FALSE))
   }
-  results <- removeDuplicateVarRanges(results)
-
-  # Remove RHSonly by passing through declRules.
-  results <- flatten(lapply(results, \(vr)
-  lapply(modelDef$declRules[[getVarName(vr)]]$rules, \(rule) {
-    nodeRange <- rule$apply(vr)
-    if (is.null(nodeRange)) {
-      return(NULL)
-    } else {
-      return(nodeRange$toVarRange(fromStochRule = vr$fromStochRule))
-    }
-  })))
+  
   if (!length(results)) {
     return(NULL)
   }
 
   if (.sort) {
+    # TODO: what happens here if have RHSonly?
     # Ordering is only relevant at calcRange stage and a single nodeRange can contain
     # elements with various sortIDs, so we convert to nodeChars first and then get their
     # sortID by creating a temporary calcRange for each.
-    nodeChars <- unlist(lapply(results, \(vr)
-    lapply(
-      modelDef$calcRules[[getVarName(vr)]]$rules,
-      \(rule) {
-        tmp <- rule$apply(vr)
-        if (!is.null(tmp)) {
-          return(tmp$toNodeChars())
-        } else {
-          return(NULL)
-        }
-      }
-    )))
-    calcRanges <- flatten(lapply(nodeChars, function(node) {
-      lapply(modelDef$calcRules[[getVarName(node)]]$rules, function(rule) {
+    nodeChars <- unlist(lapply(results, \(nr) nr$toNodeChars()))
+    calcRanges <- lapply(nodeChars, function(node) {
+      calcRange <- flatten(lapply(model$modelDef$calcRules[[getVarName(node)]]$rules, function(rule) {
         rule$makeCalcRange(rule$apply(node))
-      })
-    }))
+      }))
+      if(length(calcRange) > 1) stop("unexpected multiple calcRange results from single nodeRange")
+      if(length(calcRange) == 1) calcRange <- calcRange[[1]]
+      return(calcRange)
+    })
     if (length(nodeChars) != length(calcRanges)) {
       stop("unexpected mismatch between node character representation and calcRanges in `getNodes` sorting")
     }
-    ord <- order(sapply(calcRanges, \(x) x$sortID))
+    ord <- order(sapply(calcRanges, \(x) {
+      if(inherits(x,'calcRangeClass')) {
+        return(x$sortID)
+      } else {
+        return(-Inf) # This handles RHSonly.
+      }
+    }))
     results <- nodeChars[ord]
-    names(results) <- NULL
     if (returnScalarComponents) {
       results <- unlist(lapply(results, \(x) varRangeClass$new(x)$toVarChars(expandScalars = TRUE)))
     }
+    names(results) <- NULL
   } else {
     if (nodesAsChars) {
-      return(unlist(lapply(results, \(x) x$toVarChars(expandScalars = returnScalarComponents))))
+      if(returnScalarComponents) {
+        return(unlist(lapply(results, \(x) x$toVarChars(expandScalars = TRUE))))
+      } else {
+        return(unlist(lapply(results, \(x) x$toNodeChars())))
+      }
     } else {
       if (returnScalarComponents) { # TODO: put into new messaging system
         warning("one must request result as characters via `nodesAsChars` in order to use `returnScalarComponents`")
diff --git a/nimbleModel/tests/testthat/test-MCMCconfiguration.R b/nimbleModel/tests/testthat/test-MCMCconfiguration.R
index 4645f94..a5d85eb 100644
--- a/nimbleModel/tests/testthat/test-MCMCconfiguration.R
+++ b/nimbleModel/tests/testthat/test-MCMCconfiguration.R
@@ -2,7 +2,6 @@
 ## how to create and add to an MCMC configuration, using
 ## a nodeRange or varRange.
 
-library(nimbleModel)
 code <- quote({
     for(i in 1:5) {
         y[i] ~ dnorm(mu, sd = sigma)
diff --git a/nimbleModel/tests/testthat/test-modelGraph.R b/nimbleModel/tests/testthat/test-modelGraph.R
index 036dd09..2b9f947 100644
--- a/nimbleModel/tests/testthat/test-modelGraph.R
+++ b/nimbleModel/tests/testthat/test-modelGraph.R
@@ -1075,10 +1075,9 @@ test_that("nested indexing", {
         for(i in 1:10) 
             mu[i] ~ dnorm(0, 1)
     })
-    modelDef <- modelDefClass$new(code, constants = list(k = c(2, 4, 5), block = c(3, 5, 7, 6, 6)))
-    result <- getDependencies(modelDef, 'mu[6]', self = FALSE)
-    expect_equal(result[[1]], varRangeClass$new(list(newIndexRange(quote(2:3))), varName = 'y',
-                                                fromStochRule = TRUE))
+    model <- nimbleModel(code, constants = list(k = c(2, 4, 5), block = c(3, 5, 7, 6, 6)))
+    result <- model$getDependencies('mu[6]', self = FALSE, nodesAsChars = TRUE)
+    expect_equal(result, c('y[2]', 'y[3]'))
 
     ## Internal index constant, external index dynamic.
     code <- quote({
@@ -1089,10 +1088,9 @@ test_that("nested indexing", {
         for(i in 1:5)
             block[i] ~ dcat(p[1:10])
     })
-    modelDef <- modelDefClass$new(code, constants = list(k = c(2, 4, 5)))
-    result <- getDependencies(modelDef, 'mu[6]', self = FALSE)
-    expect_equal(result[[1]], varRangeClass$new(list(newIndexRange(quote(1:3))), varName = 'y',
-                                                fromStochRule = TRUE))
+    model <- nimbleModel(code, constants = list(k = c(2, 4, 5)))
+    result <- model$getDependencies('mu[6]', self = FALSE, nodesAsChars = TRUE)
+    expect_equal(result, paste0("y[", 1:3, "]"))
 
     ## External index constant, internal index dynamic.
     code <- quote({
@@ -1103,7 +1101,7 @@ test_that("nested indexing", {
         for(i in 1:3)
             k[i] ~ dcat(p[1:5])
     })
-    expect_error(modelDef <- modelDefClass$new(code, constants = list(block = c(3, 5, 7, 6, 6))),
+    expect_error(model <- nimbleModel(code, constants = list(block = c(3, 5, 7, 6, 6))),
                  "dynamic indexing of constants is not allowed")
 
     ##  Both internal and external indices dynamic.
@@ -1118,7 +1116,8 @@ test_that("nested indexing", {
         for(i in 1:5)
             block[i] ~ dcat(p[1:10])
     })
-    modelDef <- modelDefClass$new(code)
+    model <- nimbleModel(code)
+    modelDef <- model$modelDef
 
     expect_identical(length(modelDef$declInfo[[1]]$dynamicIndexInfo), 2L)
     expect_false(modelDef$varInfo$y$anyDynamicallyIndexed)
@@ -1126,9 +1125,8 @@ test_that("nested indexing", {
     expect_true(modelDef$varInfo$mu$anyDynamicallyIndexed)
     expect_true(modelDef$varInfo$block$anyDynamicallyIndexed)
     
-    result <- getDependencies(modelDef, 'mu[6]', self = FALSE)
-    expect_equal(result[[1]], varRangeClass$new(list(newIndexRange(quote(1:3))), varName = 'y',
-                                                fromStochRule = TRUE))
+    result <- model$getDependencies('mu[6]', self = FALSE, nodesAsChars = TRUE)
+    expect_equal(result, paste0("y[", 1:3, "]"))
 
     code <- quote({
         for(i in 1:3) 
@@ -1140,12 +1138,11 @@ test_that("nested indexing", {
         for(i in 1:5)
             block[i] ~ dcat(p[1:10])
     })
-    modelDef <- modelDefClass$new(code, data = list(block = c(3, 5, 6, 5, 2)))
+    model <- nimbleModel(code, data = list(block = c(3, 5, 6, 5, 2)))
 
     ## If `block` is not constant, we still have full set of dependencies set up.
-    result <- getDependencies(modelDef, 'mu[9]', self = FALSE)
-    expect_equal(result[[1]], varRangeClass$new(list(newIndexRange(quote(1:3))), varName = 'y',
-                                                fromStochRule = TRUE))
+    result <- model$getDependencies('mu[9]', self = FALSE, nodesAsChars = TRUE)
+    expect_equal(result, paste0('y[', 1:3 , ']'))
 
     code <- quote({
         for(i in 1:3) 
@@ -1153,11 +1150,10 @@ test_that("nested indexing", {
         for(i in 1:10) 
             mu[i] ~ dnorm(0, 1)
     })
-    expect_message(modelDef <- modelDefClass$new(code, dimensions = list(block=5)),
+    expect_message(model <- nimbleModel(code, dimensions = list(block=5)),
                                                  "Detected use of non-constant")
-    result <- getDependencies(modelDef, 'mu[6]', self = FALSE)
-    expect_equal(result[[1]], varRangeClass$new(list(newIndexRange(quote(1:3))), varName = 'y',
-                                                fromStochRule = TRUE))
+    result <- model$getDependencies('mu[6]', self = FALSE, nodesAsChars = TRUE)
+    expect_equal(result, paste0('y[', 1:3 , ']'))
                          
 })
 
@@ -1204,30 +1200,30 @@ test_that("basic check of graph interface", {
                      c('theta'))
     result <- getNodes(model, c('theta','y'))
     expect_identical(sapply(result, function(node) node$varName),
-                     c('y','theta'))
+                     c('theta','y'))
 
-    result <- getDependencies(model$modelDef, 'mu0')
+    result <- model$getDependencies('mu0')
     expect_identical(sapply(result, function(node) node$varName),
                      c('mu0','theta','y'))
-    result <- getDependencies(model$modelDef, 'mu0', self = FALSE)
+    result <- model$getDependencies('mu0', self = FALSE)
     expect_identical(sapply(result, function(node) node$varName),
                      c('theta','y'))
     
-    result <- getDependencies(model$modelDef, 'mu0', immediateOnly = TRUE)
+    result <- model$getDependencies('mu0', immediateOnly = TRUE)
     expect_identical(sapply(result, function(node) node$varName),
                      c('mu0','theta'))
 
-    result <- getDependencies(model$modelDef, c('mu0','theta'), self = FALSE)
+    result <- model$getDependencies(c('mu0','theta'), self = FALSE)
     expect_identical(sapply(result, function(node) node$varName),
                      c('y'))
 
-    result <- getParents(model$modelDef, 'y')
+    result <- model$getParents('y')
     expect_identical(sapply(result, function(node) node$varName),
                      c('theta','mu0'))
-    result <- getParents(model$modelDef, 'y', self = TRUE)
+    result <- model$getParents('y', self = TRUE)
     expect_identical(sapply(result, function(node) node$varName),
                      c('y','theta','mu0'))
-    result <- getParents(model$modelDef, 'y', immediateOnly = TRUE)
+    result <- model$getParents('y', immediateOnly = TRUE)
     expect_identical(sapply(result, function(node) node$varName),
                      c('theta'))
 
@@ -1241,21 +1237,21 @@ test_that("basic check of graph interface", {
         theta ~ dnorm(mu0, 1)
         mu0 ~ dnorm(0, sd = sigma)        
     })
-    modelDef <- modelDefClass$new(code)
+    model <- nimbleModel(code)
         
-    result <- getDependencies(modelDef, 'mu0')
+    result <- model$getDependencies('mu0')
     expect_identical(sapply(result, function(node) node$varName),
                      c('mu0','theta'))
 
-    result <- getDependencies(modelDef, 'mu0', downstream = TRUE)
+    result <- model$getDependencies('mu0', downstream = TRUE)
     expect_identical(sapply(result, function(node) node$varName),
                      c('mu0','theta','y'))
 
-    result <- getParents(modelDef, 'y')
+    result <- model$getParents('y')
     expect_identical(sapply(result, function(node) node$varName),
                      c('theta'))
 
-    result <- getParents(modelDef, 'y', upstream = TRUE)
+    result <- model$getParents('y', upstream = TRUE)
     expect_identical(sapply(result, function(node) node$varName),
                      c('theta','mu0'))
 
@@ -1265,9 +1261,9 @@ test_that("getDependencies deals with repeated children", {
     code <- quote({
         y ~ dnorm(mu, sd = tau)
     })
-    modelDef <- modelDefClass$new(code)
+    model <- nimbleModel(code)
     
-    result <- getDependencies(modelDef, c('tau','mu'))
+    result <- model$getDependencies(c('tau','mu'))
     expect_identical(sapply(result, function(node) node$varName),
                      c('y'))
 
@@ -1280,12 +1276,13 @@ test_that("getDependencies deals with repeated children", {
         tau ~ dunif(0, 1)
         sigma ~ dunif(0, 1)
     })
-    modelDef <- modelDefClass$new(code)
+    model <- nimbleModel(code)
     
-    getDependencies(modelDef, c('mu0','sigma','tau'), self = FALSE)
+    expect_identical(model$getDependencies(c('mu0','sigma','tau'), self = FALSE, nodesAsChars = TRUE),
+                     c(paste0('mu[', 1:5, ']'), paste0('y[', 1:5, ']')))
 })
 
-test_that("getDependencies deals with repeated parents", {
+test_that("getParents deals with repeated parents", {
     code <- quote({
         y ~ dnorm(mu, sd = tau)
         z ~ dnorm(mu, 1)
@@ -1293,9 +1290,9 @@ test_that("getDependencies deals with repeated parents", {
         tau ~ dunif(0,1)
         mu ~ dunif(0,1)
     })
-    modelDef <- modelDefClass$new(code)
+    model <- nimbleModel(code)
     
-    result <- getParents(modelDef, c('y','z'))
+    result <- model$getParents(c('y','z'))
     expect_identical(sapply(result, function(node) node$varName),
                      c('mu','tau'))
 
@@ -1308,13 +1305,13 @@ test_that("getDependencies with multiple children", {
         z ~ dnorm(mu, 1)
         w <- z + 3
     })
-    modelDef <- modelDefClass$new(code)
+    model <- nimbleModel(code)
 
-    result <- getDependencies(modelDef, 'mu')
+    result <- model$getDependencies('mu')
     expect_identical(sapply(result, function(node) node$varName),
                      c('y','z'))
     
-    result <- getDependencies(modelDef, 'mu', downstream = TRUE)
+    result <- model$getDependencies('mu', downstream = TRUE)
     expect_identical(sapply(result, function(node) node$varName),
                      c('y','z','w'))
 
@@ -1329,9 +1326,9 @@ test_that("getParents traversal with mix of stoch/determ edges", {
         x <- phi + 2
         phi ~ dunif(0,1)
     })
-    modelDef <- modelDefClass$new(code)
+    model <- nimbleModel(code)
     
-    result <- getParents(modelDef, 'y')
+    result <- model$getParents('y')
     expect_identical(sapply(result, function(node) node$varName),
                      c('theta','z','w','x','phi'))
 
@@ -1342,9 +1339,9 @@ test_that("same dependent on RHS", {
         y ~ dnorm(mu, sd = mu)
         mu ~ dunif(0,1)
     })
-    modelDef <- modelDefClass$new(code)
+    model <- nimbleModel(code)
 
-    result <- getParents(modelDef, 'y')
+    result <- model$getParents('y')
     expect_identical(sapply(result, function(node) node$varName),
                     'mu')
 })
@@ -1375,7 +1372,7 @@ test_that("basic hierarchical models", {
 
     result <- getNodes(model, c('z[1:5]', 'mu'))
     expect_identical(sapply(result, function(node) node$varName),
-                     c('mu','z'))
+                     c('z','mu'))
 
     result <- getNodes(model, topOnly = TRUE)
     expect_identical(sapply(result, function(node) node$varName),
@@ -1404,45 +1401,45 @@ test_that("basic hierarchical models", {
     expect_identical(sapply(result, function(node) node$varName),
                      c('y','mu','tau','sigma','z','mu0','bnd'))
 
-    result <- getDependencies(model$modelDef, 'sigma')
+    result <- getDependencies(model, 'sigma')
     expect_identical(sapply(result, function(node) node$varName),
                      c('sigma','mu'))
     expect_equal(result[[2]]$indexRanges, 
                list(newIndexRange(quote(1:10))))
 
-    result <- getDependencies(model$modelDef, 'mu', self = FALSE)
+    result <- getDependencies(model, 'mu', self = FALSE)
     expect_identical(sapply(result, function(node) node$varName),
                      c('y'))
     expect_equal(result[[1]]$indexRanges, 
                list(newIndexRange(quote(1:10))))
     
-    result <- getDependencies(model$modelDef, 'y', self = FALSE)
+    result <- getDependencies(model, 'y', self = FALSE)
     expect_equal(result[[1]]$indexRanges, 
                list(newIndexRange(matrix(k))))
     
-    result <- getDependencies(model$modelDef, 'y[1:3]', self = FALSE)
+    result <- getDependencies(model, 'y[1:3]', self = FALSE)
     expect_equal(result[[1]]$indexRanges, 
                list(newIndexRange(2)))
 
-    result <- getDependencies(model$modelDef, 'sigma', downstream = TRUE)
+    result <- getDependencies(model, 'sigma', downstream = TRUE)
     expect_identical(sapply(result, function(node) node$varName),
                      c('sigma','mu','y','z'))
 
-    result <- getParents(model$modelDef, 'y')
+    result <- getParents(model, 'y')
     expect_identical(sapply(result, function(node) node$varName),
                      c('mu','tau'))
 
-    result <- getParents(model$modelDef, 'y', upstream = TRUE)
+    result <- getParents(model, 'y', upstream = TRUE)
     expect_identical(sapply(result, function(node) node$varName),
                      c('mu','tau','sigma'))
 
-    result <- getParents(model$modelDef, 'y[1:5]', upstream = TRUE)
+    result <- getParents(model, 'y[1:5]', upstream = TRUE)
     expect_identical(sapply(result, function(node) node$varName),
                      c('mu','tau','sigma'))
     expect_equal(result[[1]]$indexRanges, 
                list(newIndexRange(quote(1:5))))
 
-    result <- getParents(model$modelDef, varRangeClass$new(list(newIndexRange(quote(1:5))), varName = 'y'),
+    result <- getParents(model, varRangeClass$new(list(newIndexRange(quote(1:5))), varName = 'y'),
                          upstream = TRUE)
     expect_identical(sapply(result, function(node) node$varName),
                      c('mu','tau','sigma'))
@@ -1477,7 +1474,7 @@ test_that("basic hierarchical models", {
     expect_equal(result[[1]]$indexRanges, 
                list(newIndexRange(quote(1:10)), newIndexRange(quote(1:10))))
 
-    result <- getDependencies(model$modelDef, 'pr[2,2]', self = FALSE)
+    result <- getDependencies(model, 'pr[2,2]', self = FALSE)
     expect_identical(sapply(result, function(node) node$varName),
                      c('w','mu'))
     expect_equal(result[[1]]$indexRanges, 
@@ -1485,14 +1482,14 @@ test_that("basic hierarchical models", {
     expect_equal(result[[2]]$indexRanges, 
                list(newIndexRange(quote(1:10))))
 
-    result <- getDependencies(model$modelDef, 'pr[2,2]', self = FALSE, downstream = TRUE)
+    result <- getDependencies(model, 'pr[2,2]', self = FALSE, downstream = TRUE)
     expect_identical(sapply(result, function(node) node$varName),
                      c('w','mu','y'))
 
-    result <- getParents(model$modelDef, 'pr[2,2]')
+    result <- getParents(model, 'pr[2,2]')
     expect_identical(result, NULL)
 
-    result <- getParents(model$modelDef, 'mu[3]')
+    result <- getParents(model, 'mu[3]')
     expect_identical(sapply(result, function(node) node$varName),
                      c('mu0','pr','mu00'))
     expect_equal(result[[2]]$indexRanges, 
@@ -1523,7 +1520,7 @@ test_that("basic hierarchical models", {
     expect_equal(result[[1]]$indexRanges, 
                list(newIndexRange(quote(2))))
 
-    result <- getDependencies(model$modelDef, 'beta[2]')
+    result <- getDependencies(model, 'beta[2]')
     expect_identical(sapply(result, function(node) node$varName),
                      c('beta','mn','y'))
     expect_equal(result[[1]]$indexRanges, 
@@ -1533,15 +1530,15 @@ test_that("basic hierarchical models", {
     expect_equal(result[[3]]$indexRanges, 
                list(newIndexRange(quote(1:5)), newIndexRange(quote(1:3))))
 
-    result <- getParents(model$modelDef, 'y[1,1:3]')
+    result <- getParents(model, 'y[1,1:3]')
     expect_identical(sapply(result, function(node) node$varName),
                      c('mn','pr','beta'))
 
-    result <- getParents(model$modelDef, 'y[1,2]')
+    result <- getParents(model, 'y[1,2]')
     expect_identical(sapply(result, function(node) node$varName),
                      c('mn','pr','beta'))
     
-    result <- getParents(model$modelDef, c('y[1,2]','y[2,3]'))
+    result <- getParents(model, c('y[1,2]','y[2,3]'))
     expect_identical(sapply(result, function(node) node$varName),
                      c('mn','pr','beta'))
     
@@ -1553,8 +1550,8 @@ test_that("mixed-length block dependences", {
         for(i in 1:3)
             y[i, n1[i]:n2[i]] ~ dmulti(p[n1[i]:n2[i]], 10)
     })
-    modelDef <- modelDefClass$new(code, constants = list(n1 = c(3,1,2), n2 = c(6,3,2)))
-    result <- getDependencies(modelDef, 'p[2]')
+    model <- nimbleModel(code, constants = list(n1 = c(3,1,2), n2 = c(6,3,2)))
+    result <- getDependencies(model, 'p[2]')
     expect_equal(result[[1]]$indexRanges, list(newIndexRange(matrix(c(2,2,2,3,1,2,3,2), ncol = 2))))
 })
 
@@ -1612,28 +1609,17 @@ test_that("state-space model", {
     expect_equal(result[[1]]$indexRanges, 
                list(newIndexRange(3)))
 
-    result <- getDependencies(modelDef, 'y[1]')
-    expect_length(result, 1)
-    expect_equal(result[[1]],
-                 varRangeClass$new(list(newIndexRange(2)), varName = 'y', fromStochRule = TRUE))
+    result <- getDependencies(model, 'y[1]', nodesAsChars = TRUE)
+    expect_identical(result, 'y[2]')
     
-    result <- getDependencies(modelDef, 'y[2]')
-    expect_length(result, 2)
+    result <- getDependencies(model, 'y[2]', nodesAsChars = TRUE)
+    expect_identical(result, c('y[2]', 'y[3]'))
 
-    expect_equal(result[[1]],
-                 varRangeClass$new(list(newIndexRange(2)), varName = 'y', fromStochRule = TRUE))
-    expect_equal(result[[2]],
-                 varRangeClass$new(list(newIndexRange(3)), varName = 'y'))
+    result <- getDependencies(model, 'y[2]', downstream = TRUE, nodesAsChars = TRUE)
+    expect_identical(result, c('y[2]', 'y[3]', 'y[4]'))
 
-    result <- getDependencies(modelDef, 'y[2]', downstream = TRUE)
-    expect_length(result, 3)
-    expect_equal(result[[3]],
-                 varRangeClass$new(list(newIndexRange(4)), varName = 'y'))
-
-    result <- getParents(modelDef, 'y[3]')
-    expect_length(result, 1)
-    expect_equal(result[[1]],
-                 varRangeClass$new(list(newIndexRange(2)), varName = 'y'))
+    result <- getParents(model, 'y[3]', nodesAsChars = TRUE)
+    expect_identical(result, 'y[2]')
 
     code <- quote({
         for(i in 2:4)
@@ -1641,12 +1627,9 @@ test_that("state-space model", {
         y[1] ~ dnorm(0,1)
     })
     model <- nimbleModel(code)
-    modelDef <- model$modelDef
     
-    result <- getParents(modelDef, 'y[3]', upstream = TRUE)
-    expect_length(result, 2)
-    expect_equal(result[[2]],
-                 varRangeClass$new(list(newIndexRange(1)), varName = 'y'))
+    result <- getParents(model, 'y[3]', upstream = TRUE, nodesAsChars = TRUE)
+    expect_identical(result, c('y[2]','y[1]'))
 
     code <- quote({
         for(i in 1:5)
@@ -1658,7 +1641,6 @@ test_that("state-space model", {
         sigma ~ dunif(0, 1)
     })
     model <- nimbleModel(code)
-    modelDef <- model$modelDef
 
     result <- getNodes(model, topOnly = TRUE)
     expect_length(result, 3)
@@ -1671,39 +1653,29 @@ test_that("state-space model", {
     expect_length(result, 1)
     expect_identical(result[[1]]$varName, 'z')
 
-    result <- getDependencies(modelDef, 'z')
+    result <- getDependencies(model, 'z')
     expect_length(result, 3)
     expect_identical(sapply(result, function(node) node$varName),
                      c('z','z','y'))
-    expect_equal(result[[1]]$indexRanges, 
-                     list(newIndexRange(quote(2:5))))
-    expect_equal(result[[2]]$indexRanges, 
-                     list(newIndexRange(1)))
+    expect_identical(result[[1]]$toNodeChars(), paste0('z[', 2:5, ']'))
+    expect_identical(result[[2]]$toNodeChars(), 'z[1]')
+    expect_identical(result[[3]]$toNodeChars(), paste0('y[', 1:5, ']'))
 
-    result <- getDependencies(modelDef, 'z', self = FALSE)
+    result <- getDependencies(model, 'z', self = FALSE)
+    expect_identical(result[[1]]$toNodeChars(), paste0('y[', 1:5, ']'))
     expect_length(result, 1)
-    expect_identical(sapply(result, function(node) node$varName),
-                     c('y'))
     
-    result <- getDependencies(modelDef, 'z[3]')
-    expect_length(result, 3)
-    expect_identical(sapply(result, function(node) node$varName),
-                     c('z','z','y'))
-    expect_equal(result[[1]]$indexRanges, 
-                     list(newIndexRange(quote(3))))
-    expect_equal(result[[2]]$indexRanges, 
-                     list(newIndexRange(quote(4))))
-    expect_equal(result[[3]]$indexRanges, 
-                     list(newIndexRange(quote(3))))
+    result <- getDependencies(model, 'z[3]')
+    expect_length(result, 2)
+    expect_identical(result[[1]]$toNodeChars(), paste0('z[', 3:4, ']'))
+    expect_identical(result[[2]]$toNodeChars(), 'y[3]')
+    
 
-    result <- getDependencies(modelDef, 'z[3]', self = FALSE)
+    result <- getDependencies(model, 'z[3]', self = FALSE)
     expect_length(result, 2)
-    expect_identical(sapply(result, function(node) node$varName),
-                     c('y','z'))
-    expect_equal(result[[1]]$indexRanges, 
-                     list(newIndexRange(quote(3))))
-    expect_equal(result[[2]]$indexRanges, 
-                     list(newIndexRange(quote(4))))
+    expect_identical(result[[1]]$toNodeChars(), 'y[3]')
+    expect_identical(result[[2]]$toNodeChars(), 'z[4]')
+
 
 })
 
@@ -1714,15 +1686,14 @@ test_that("error trapping for unexpected vars/nodes", {
         y[i]~dnorm(y[i-1],1)
     })
     model <- nimbleModel(code)
-    modelDef <- model$modelDef
 
     expect_null(getNodes(model, 'x'))
     expect_null(getNodes(model, 'x[3]'))
     expect_null(getNodes(model, 'y[20]'))
     
-    expect_null(getDependencies(modelDef, 'x'))
-    expect_null(getDependencies(modelDef, 'x[1]'))
-    expect_null(getDependencies(modelDef, 'y[20]'))
+    expect_null(getDependencies(model, 'x'))
+    expect_null(getDependencies(model, 'x[1]'))
+    expect_null(getDependencies(model, 'y[20]'))
 
 })
 
@@ -1763,7 +1734,6 @@ test_that("complicated input varRange", {
                     y[i+1,j,k] ~ dnorm(theta[i,j,k], 1)
     })
     model <- nimbleModel(code)
-    modelDef <- model$modelDef
 
     vr <- varRangeClass$new(list(newIndexRange(quote(2:3)),
                                  newIndexRange(matrix(c(2,4,5,1), nrow = 2))),
@@ -1778,15 +1748,10 @@ test_that("complicated input varRange", {
     expect_equal(result[[1]]$indexRanges,
                      vr$indexRanges[c(2,1)])
     
-    result <- getDependencies(modelDef, vr)
+    result <- getDependencies(model, vr)
     expect_length(result, 1)
-    expect_identical(sapply(result, function(node) node$varName),
-                     c('y'))
-    expect_equal(result[[1]],
-                 varRangeClass$new(list(newIndexRange(matrix(c(3,5,5,1), nrow = 2)),
-                                        newIndexRange(quote(2:3))),
-                            rangeToIndexSlot = list(c(1,3), 2),
-                            varName = 'y', fromStochRule = TRUE))
+    expect_identical(result[[1]]$toNodeChars(),
+                     c("y[5, 2, 1]", "y[3, 2, 5]", "y[5, 3, 1]", "y[3, 3, 5]"))
 })
 
 test_that("duplicated RHS elements", {
@@ -1795,18 +1760,13 @@ test_that("duplicated RHS elements", {
         mu ~ dnorm(0, 1)
     })
     
-    modelDef <- modelDefClass$new(code)
-    expect_equal(getParents(modelDef, 'y')[[1]],
-                 varRangeClass$new('mu', fromStochRule = TRUE))
+    model <- nimbleModel(code)
+    expect_identical(getParents(model, 'y', nodesAsChars = TRUE), "mu")
 
-    deps <- getDependencies(modelDef, 'mu')
+    deps <- getDependencies(model, 'mu')
     expect_identical(length(deps), 2L)
-    expect_equal(deps[[1]],
-                 varRangeClass$new('mu'))
-    expect_equal(deps[[2]],
-                 varRangeClass$new('y', fromStochRule = TRUE))
-    
-
+    expect_equal(deps[[1]]$toNodeChars(), 'mu')
+    expect_equal(deps[[2]]$toNodeChars(), 'y')
 })
 
 test_that("handling of `self` in graph traversal", {
@@ -1819,25 +1779,28 @@ test_that("handling of `self` in graph traversal", {
         theta ~ dnorm(0, 1)
     })
     
-    modelDef <- modelDefClass$new(code)
-
-    deps <- getDependencies(modelDef, c('mu[1:5]','theta'), self = TRUE)
-    expect_equal(deps, list(varRangeClass$new(list(), varName = 'theta'),
-                            varRangeClass$new(list(newIndexRange(quote(1:5))), varName = 'mu', fromStochRule = TRUE),
-                            varRangeClass$new(list(newIndexRange(quote(1:5))), varName = 'y', fromStochRule = TRUE)))
+    model <- nimbleModel(code)
 
-    deps <- getDependencies(modelDef, c('mu[1:5]','theta'), self = FALSE)
+    deps <- getDependencies(model, c('mu[1:5]','theta'), self = TRUE)
+    expect_length(deps, 3)
+    expect_identical(deps[[1]]$toNodeChars(), paste0('mu[', 1:5, ']'))
+    expect_identical(deps[[2]]$toNodeChars(), 'theta')
+    expect_identical(deps[[3]]$toNodeChars(), paste0('y[', 1:5, ']'))
+ 
+    deps <- getDependencies(model, c('mu[1:5]','theta'), self = FALSE)
+    expect_length(deps, 1)
+    expect_identical(deps[[1]]$toNodeChars(), paste0('y[', 1:5, ']'))
 
-    expect_equal(deps, list(varRangeClass$new(list(newIndexRange(quote(1:5))), varName = 'y', fromStochRule = TRUE)))
-    deps <- getDependencies(modelDef, c('mu[1:3]','theta'), self = FALSE)
-    expect_equal(deps, list(varRangeClass$new(list(newIndexRange(quote(1:3))), varName = 'y', fromStochRule = TRUE),
-                            varRangeClass$new(list(newIndexRange(quote(4:5))), varName = 'mu')))
+    deps <- getDependencies(model, c('mu[1:3]','theta'), self = FALSE)
+    expect_length(deps, 2)
+    expect_identical(deps[[1]]$toNodeChars(), paste0('y[', 1:3, ']'))
+    expect_identical(deps[[2]]$toNodeChars(), paste0('mu[', 4:5, ']'))
     
-    deps <- getDependencies(modelDef, c('mu[1:3]','theta'), self = TRUE)
-    expect_equal(deps, list(varRangeClass$new(list(), varName = 'theta'),
-                            varRangeClass$new(list(newIndexRange(quote(1:3))), varName = 'mu', fromStochRule = TRUE),
-                            varRangeClass$new(list(newIndexRange(quote(4:5))), varName = 'mu'),
-                            varRangeClass$new(list(newIndexRange(quote(1:3))), varName = 'y', fromStochRule = TRUE)))
+    deps <- getDependencies(model, c('mu[1:3]','theta'), self = TRUE)
+    expect_length(deps, 3)
+    expect_identical(deps[[1]]$toNodeChars(), paste0('mu[', 1:5, ']'))
+    expect_identical(deps[[2]]$toNodeChars(), 'theta')
+    expect_identical(deps[[3]]$toNodeChars(), paste0('y[', 1:3, ']'))
 
 })
 
@@ -2525,17 +2488,16 @@ test_that("SSM with additional pieces handled", {
         z2[1] ~ dnorm(y2[1], 1)
         
     })
-    modelDef <- modelDefClass$new(code)
-    
+    model <- nimbleModel(code)
     
-    expect_equal(getDependencies(modelDef, 'sigma', self = FALSE),
-                 list(varRangeClass$new(list(newIndexRange(quote(1:5))), varName = 'y',
-                                        fromStochRule = TRUE),
-                      varRangeClass$new(list(newIndexRange(quote(1:5))), varName = 'y2',
-                                        fromStochRule = TRUE)))
-    expect_equal(modelDef$endRules$y2$rules[[1]]$fullRange,
+    result <- getDependencies(model, 'sigma', self = FALSE)
+    expect_length(result, 2)
+    expect_identical(result[[1]]$toNodeChars(), paste0('y[', 1:5, ']'))
+    expect_identical(result[[2]]$toNodeChars(), paste0('y2[', 1:5, ']'))
+
+    expect_equal(model$modelDef$endRules$y2$rules[[1]]$fullRange,
                  varRangeClass$new(list(newIndexRange(quote(2:5))), varName = 'y2'))
-    expect_equal(modelDef$latentRules$y2$rules[[1]]$fullRange,
+    expect_equal(model$modelDef$latentRules$y2$rules[[1]]$fullRange,
                  varRangeClass$new(list(newIndexRange(quote(1))), varName = 'y2'))
     
 })
@@ -2548,30 +2510,27 @@ test_that("duplicated nodes", {
             mu[i] ~ dnorm(0, 1)
         }
     })
-    modelDef <- modelDefClass$new(code, constants = list(k = c(1,1,2)))
-    expect_equal(getParents(modelDef, 'y')[[1]],
-                 varRangeClass$new(list(newIndexRange(quote(1:2))),
-                                   varName = 'mu', fromStochRule = TRUE))
-    expect_equal(getDependencies(modelDef, 'mu[1]', self = FALSE)[[1]],
-                 varRangeClass$new(list(newIndexRange(quote(1:2))),
-                                   varName = 'y', fromStochRule = TRUE))
+    model <- nimbleModel(code, constants = list(k = c(1,1,2)))
+    expect_identical(getParents(model, 'y', nodesAsChars = TRUE),
+                 c('mu[1]','mu[2]'))
+    expect_equal(getDependencies(model, 'mu[1]', self = FALSE, nodesAsChars = TRUE),
+                 c('y[1]','y[2]'))
 })
 
 test_that("missing indexing", {
     modelCode <- quote({
         y <- sum(mu[])
     })
-    m <- nimbleModel(modelCode, dimensions = list(mu = 5))
-    modelDef <- m$modelDef
+    model <- nimbleModel(modelCode, dimensions = list(mu = 5))
 
-    expect_equal(getNodes(m, includeRHSonly = TRUE)[[2]]$toVarRange(),
+    expect_equal(getNodes(model, includeRHSonly = TRUE)[[2]]$toVarRange(),
                  varRangeClass$new(list(newIndexRange(quote(1:5))),
                                         varName = 'mu'))
 
-    expect_equal(getDependencies(modelDef, 'mu[2]')[[1]],
-                 varRangeClass$new(list(), varName = 'y', fromStochRule = FALSE))
+    expect_identical(getDependencies(model, 'mu[2]', nodesAsChars = TRUE),
+                     "y")
     
-    expect_identical(getDependencies(modelDef, 'mu[9]'), NULL)
+    expect_identical(getDependencies(model, 'mu[9]'), NULL)
 
 })
 
@@ -2601,3 +2560,48 @@ test_that("flexible use of indexing in models", {
                  "unable to process")
     nimbleModel:::nimbleModelOptions(verbose = TRUE)
 })
+
+test_that("handling duplication and full nodes with getDeps/getParents", {
+    code <- nimbleCode({
+        for(i in 1:4) y[i] ~ dnorm(mu, 1)
+        for(i in 3:4) z[i] ~ dnorm(y[i],1)
+    })
+    m <- nimbleModel(code)
+    expect_identical(m$getDependencies(c('y[3]','y'), nodesAsChars = TRUE),
+                     c(paste0('y[', 1:4, ']'), 'z[3]', 'z[4]'))
+
+    code <- nimbleCode({
+        for(i in 1:4) y[i] ~ dnorm(mu[i],1)
+        mu[1:4] ~ dmnorm(z[1:4],pr[1:4,1:4])
+    })
+    m <- nimbleModel(code)
+
+    expect_identical(m$getDependencies('mu[1]', nodesAsChars = TRUE),
+                     c("mu[1:4]", "y[1]"))
+    
+    expect_identical(m$getParents('y[2]', nodesAsChars = TRUE), "mu[1:4]")
+    
+})
+
+test_that("handling RHSonly with getParents", {
+    code <- nimbleCode({
+        y ~ dnorm(mu, sd = sigma)
+        mu ~ dnorm(0,1)
+    })
+    m <- nimbleModel(code)
+    expect_identical(m$getParents('y', nodesAsChars = TRUE), 'mu')
+    expect_identical(m$getParents('y', includeRHSonly = TRUE, nodesAsChars = TRUE),
+                     c('mu', 'sigma'))
+    
+    code <- nimbleCode({
+        for(i in 1:4) y[i] ~ dnorm(mu[i],1)
+    })
+    m <- nimbleModel(code)
+    expect_identical(m$getParents('y[2]', includeRHSonly = TRUE, nodesAsChars = TRUE), "mu[2]")
+
+    expect_identical(m$getParents('y[2]', includeRHSonly = TRUE, nodesAsChars = TRUE, .sort = TRUE),
+                     c("mu[2]"))
+    expect_identical(m$getParents('y[2]', self = TRUE, includeRHSonly = TRUE, nodesAsChars = TRUE, .sort = TRUE),
+                     c("mu[2]","y[2]"))
+    
+})
diff --git a/nimbleModel/tests/testthat/test-nodeChars.R b/nimbleModel/tests/testthat/test-nodeChars.R
index 377f434..28b1158 100644
--- a/nimbleModel/tests/testthat/test-nodeChars.R
+++ b/nimbleModel/tests/testthat/test-nodeChars.R
@@ -11,19 +11,16 @@ test_that("use of nodes as characters", {
   m <- nimbleModel(code)
   nodeRanges <- m$getNodes()
   expect_true(all(sapply(nodeRanges, \(x) inherits(x, 'nodeRangeClass'))))
+
   setNimbleModelOption('nodesAsChars', TRUE)
+
   chars <- m$getNodes()
   expect_identical(chars, c("y[1, 1]","y[2, 1]","y[1, 2]","y[2, 2]","y[1, 3]","y[2, 3]","mu"))
   chars <- m$getNodes(returnScalarComponents = TRUE)
   expect_identical(chars, c("y[1, 1]","y[2, 1]","y[1, 2]","y[2, 2]","y[1, 3]","y[2, 3]","mu"))
 
   deps <- m$getDependencies('mu')
-  expect_identical(deps, c("mu", "y[1:2, 1:3]"))
-  deps <- m$getDependencies('mu',returnScalarComponents=TRUE)
   expect_identical(deps, c("mu","y[1, 1]","y[2, 1]","y[1, 2]","y[2, 2]","y[1, 3]","y[2, 3]"))
-  setNimbleModelOption('nodesAsChars', FALSE)
-  varRanges <- m$getDependencies('mu')
-  expect_true(all(sapply(varRanges, \(x) inherits(x, 'varRangeClass'))))
 
   code <- nimbleCode({
     for(i in 1:3)
@@ -31,18 +28,14 @@ test_that("use of nodes as characters", {
     for(j in 1:2)
       mu[j] ~dnorm(0,1)
   })
+  m <- nimbleModel(code)
   
-  m <- nimbleModel(code, inits = list(prec=diag(2)))
-  nodeRanges <- m$getNodes()
-  expect_true(all(sapply(nodeRanges, \(x) inherits(x, 'nodeRangeClass'))))
-  
-  setNimbleModelOption('nodesAsChars', TRUE)
   chars <- m$getNodes()
   expect_identical(chars, c("lifted_chol_oPprec_oB1to2_comma_1to2_cB_cP[1:2, 1:2]","y[1, 1:2]","y[2, 1:2]", "y[3, 1:2]", "mu[1]", "mu[2]"))
   chars <- m$getNodes(returnScalarComponents = TRUE)
   expect_identical(chars, c("lifted_chol_oPprec_oB1to2_comma_1to2_cB_cP[1, 1]","lifted_chol_oPprec_oB1to2_comma_1to2_cB_cP[2, 1]","lifted_chol_oPprec_oB1to2_comma_1to2_cB_cP[1, 2]","lifted_chol_oPprec_oB1to2_comma_1to2_cB_cP[2, 2]","y[1, 1]","y[2, 1]","y[3, 1]","y[1, 2]","y[2, 2]","y[3, 2]","mu[1]","mu[2]"))
   deps <- m$getDependencies('mu')
-  expect_identical(deps, c("mu[1:2]", "y[1:3, 1:2]"))
+  expect_identical(deps, c("mu[1]", "mu[2]", "y[1, 1:2]", "y[2, 1:2]", "y[3, 1:2]"))
   deps <- m$getDependencies('mu',returnScalarComponents = TRUE)
   expect_identical(deps, c("mu[1]","mu[2]","y[1, 1]","y[2, 1]","y[3, 1]","y[1, 2]","y[2, 2]","y[3, 2]"))
   setNimbleModelOption('nodesAsChars', FALSE)
@@ -72,7 +65,7 @@ test_that("old model API calls", {
   expect_identical(chars, c("x", "mu", "lifted_mu_plus_x", "y[1, 1]","y[2, 1]","y[1, 2]", "y[2, 2]", "y[1, 3]", "y[2, 3]"))
 
   chars <- m$expandNodeNames(c('mu','y','x'))
-  expect_identical(chars, c("y[1, 1]","y[2, 1]","y[1, 2]", "y[2, 2]", "y[1, 3]", "y[2, 3]", "mu", "x"))
+  expect_identical(chars, c("mu", "y[1, 1]","y[2, 1]","y[1, 2]", "y[2, 2]", "y[1, 3]", "y[2, 3]", "x"))
   chars <- m$expandNodeNames(c('mu','y','x'), sort = TRUE)
   expect_identical(chars, c("x", "mu", "y[1, 1]","y[2, 1]","y[1, 2]", "y[2, 2]", "y[1, 3]", "y[2, 3]"))
   chars <- m$topologicallySortNodes(c('mu','y','x'))
@@ -190,7 +183,8 @@ test_that("old model API calls", {
   truth <- c("tau", "lifted_d1_over_sqrt_oPtau_cP", paste0("y[", 1:5, "]"))
   expect_identical(m$getDependencies(c('y','tau'), .sort = TRUE, nodesAsChars = TRUE), truth)
   expect_identical(m$getParents('y', .sort = TRUE, nodesAsChars = TRUE, self = TRUE), truth)
-  expect_identical(m$getParents('y', .sort = TRUE, nodesAsChars = TRUE, self = TRUE, returnScalarComponents = TRUE),
+  expect_identical(m$getParents('y', .sort = TRUE, nodesAsChars = TRUE, self = TRUE,
+                                returnScalarComponents = TRUE),
                    truth)
 })
 
@@ -217,9 +211,9 @@ test_that("Use of .sort in cases with multiple and/or overlapping sortID values"
     y[1] ~ dnorm(0,1)
   })
   m <- nimbleModel(code)
-  truth <- c("tau", "y[1]", "lifted_rho_times_y_oBi_minus_1_cB_L2[2]","lifted_d1_over_sqrt_oPtau_cP", "y[2]", "lifted_rho_times_y_oBi_minus_1_cB_L2[3]", "y[3]", "lifted_rho_times_y_oBi_minus_1_cB_L2[4]","y[4]" ,"lifted_rho_times_y_oBi_minus_1_cB_L2[5]","y[5]" , "lifted_rho_times_y_oBi_minus_1_cB_L2[6]","y[6]")
+  truth <- c("y[1]", "tau", "lifted_rho_times_y_oBi_minus_1_cB_L2[2]","lifted_d1_over_sqrt_oPtau_cP", "y[2]", "lifted_rho_times_y_oBi_minus_1_cB_L2[3]", "y[3]", "lifted_rho_times_y_oBi_minus_1_cB_L2[4]","y[4]" ,"lifted_rho_times_y_oBi_minus_1_cB_L2[5]","y[5]" , "lifted_rho_times_y_oBi_minus_1_cB_L2[6]","y[6]")
   expect_identical(m$getNodes(.sort=TRUE,nodesAsChars=TRUE), truth)
-  expect_identical(m$getParents('y', .sort=TRUE, nodesAsChars = TRUE, self = TRUE), truth[c(2,1,3:length(truth))])  # Not clear why getParents reverses y[1] and tau, but it's not material and will change possibly once getParents calls getNodes.
+  expect_identical(m$getParents('y', .sort=TRUE, nodesAsChars = TRUE, self = TRUE), truth)
   expect_identical(m$getParents('y[4]', .sort=TRUE, nodesAsChars = TRUE, self = TRUE), c('tau','lifted_d1_over_sqrt_oPtau_cP','y[3]','lifted_rho_times_y_oBi_minus_1_cB_L2[4]', 'y[4]'))
   nrs <- m$getNodes()
   expect_identical(m$topologicallySortNodes(nrs), truth)
diff --git a/nimbleModel/tests/testthat/test-scoping.R b/nimbleModel/tests/testthat/test-scoping.R
index 0274f3c..0655dba 100644
--- a/nimbleModel/tests/testthat/test-scoping.R
+++ b/nimbleModel/tests/testthat/test-scoping.R
@@ -61,7 +61,7 @@ test_that("avoid using indexing values from global", {
 
     p <- 3
     expect_error(getNodes(model, 'y[1:p]'), "must involve two positive")
-    expect_error(getDependencies(model$modelDef, 'y[1:p]'), "must involve two positive")
+    expect_error(getDependencies(model, 'y[1:p]'), "must involve two positive")
     
 })
    
diff --git a/nimbleModel/tests/testthat/test-setupMargNodes.R b/nimbleModel/tests/testthat/test-setupMargNodes.R
index 917455e..cdb7d62 100644
--- a/nimbleModel/tests/testthat/test-setupMargNodes.R
+++ b/nimbleModel/tests/testthat/test-setupMargNodes.R
@@ -99,11 +99,11 @@ test_that("setupMargNodes/GCIS works in model with an extra edge from param to r
   expect_identical(unlist(lapply(SMN$randomEffectsNodes, \(x) x$toNodeChars())), c("x[1]", "x[2]", "y[1]", "y[2]"))
 
   SMN <- setupMargNodes(m, randomEffectsNodes = c("y[1]", "x[2]"))
-  expect_identical(unlist(lapply(SMN$paramNodes, \(x) x$toNodeChars())), c("mu", "sigma", "x[1]"))
+  expect_identical(unlist(lapply(SMN$paramNodes, \(x) x$toNodeChars())), c("x[1]", "mu", "sigma"))
 
   SMN <- setupMargNodes(m, randomEffectsNodes = c("mu", "x[1]"))
-  expect_identical(unlist(lapply(SMN$paramNodes, \(x) x$toNodeChars())), c("sigma", "y[2]"))
-  expect_identical(unlist(lapply(SMN$givenNodes, \(x) x$toNodeChars())), c('sigma','y[2]','x[2]','y[1]','z[1]','z[2]'))
+  expect_identical(unlist(lapply(SMN$paramNodes, \(x) x$toNodeChars())), c("y[2]", "sigma"))
+  expect_identical(unlist(lapply(SMN$givenNodes, \(x) x$toNodeChars())), c('y[2]','sigma','x[2]','y[1]','z[1]','z[2]'))
 
   # Warning from missing deterministic nodes
   expect_message(SMN <- setupMargNodes(m, calcNodes = c("x[1]","x[2]","y[1]","y[2]","z[1]","z[2]")),
@@ -150,11 +150,11 @@ test_that("getConditionallyIndependentSets works in model with one set and deter
     data = list(Y1 = 1)
   )
   # All sets
-  expect_identical(getConditionallyIndependentSets(m), list(c("REA1", "REA2", "REB1", "REB2")))
+  expect_identical(getConditionallyIndependentSets(m), list(c("REA1", "REA2", "REB2", "REB1")))
   # first-stage latents stay separated
   expect_identical(getConditionallyIndependentSets(m, c("REA1", "REB1")), list("REA1", "REB1"))
   # first-stage latents are combined if unknownAsGiven is FALSE
-  expect_identical(getConditionallyIndependentSets(m, c("REA1", "REB1"), unknownAsGiven=FALSE), list(c("REA1", "REA2", "REB1", "REB2")))
+  expect_identical(getConditionallyIndependentSets(m, c("REA1", "REB1"), unknownAsGiven=FALSE), list(c("REA1", "REA2", "REB2", "REB1")))
   # second-stage latents are connected
   expect_identical(getConditionallyIndependentSets(m, c("REA2", "REB2")), list(c("REA2", "REB2")))
   # TODO: need to allow deterministic as given; currently error-trapped.
@@ -172,14 +172,14 @@ test_that("getConditionallyIndependentSets works in model with one set and deter
   SMN <- setupMargNodes(m)
   expect_identical(SMN$paramNodes[[1]]$toNodeChars(), "P1")
   expect_identical(unlist(lapply(SMN$randomEffectsNodes, \(x) x$toNodeChars())), c("REA1", "REA2", "REB1", "REB2"))
-  expect_identical(SMN$randomEffectsSets, list(c('REA1','REA2','REB1','REB2')))
+  expect_identical(SMN$randomEffectsSets, list(c('REA1','REA2','REB2','REB1')))
 
   SMN <- setupMargNodes(m, paramNodes = "REA1")
   expect_identical(SMN$randomEffectsNodes[[1]]$toNodeChars(), c('REA2'))
-  expect_identical(unlist(lapply(SMN$calcNodes, \(x) x$toNodeChars())), c('D3C3','Y1','REA2','D3'))
+  expect_identical(unlist(lapply(SMN$calcNodes, \(x) x$toNodeChars())), c('REA2','D3', 'D3C3','Y1'))
 
   SMN <- setupMargNodes(m, calcNodes = m$getDependencies("REB2"))
-  expect_identical(unlist(lapply(SMN$paramNodes, \(x) x$toNodeChars())), c("REA2", "REB1"))
+  expect_identical(unlist(lapply(SMN$paramNodes, \(x) x$toNodeChars())), c("REB1", "REA2"))
 
   SMN <- setupMargNodes(m, randomEffectsNodes = c("P1", "REA1","REA2", "REB1","REB2"))
   expect_identical(SMN$paramNodes, NULL)
@@ -206,11 +206,11 @@ test_that("getConditionallyIndependentSets works in state-space model with a cou
   expect_identical(getConditionallyIndependentSets(m, "x[1:2]", givenNodes = c("w[1]", "y"), unknownAsGiven=FALSE),
                    list(c("x[1]", "x[2]", "x[3]", "x[4]")))
   expect_identical(getConditionallyIndependentSets(m, "x[2]", givenNodes = c("w[1]", "y"), unknownAsGiven=FALSE),
-                   list(paste0("x[", 1:4, "]")))
+                   list(paste0("x[", c(2:4,1), "]")))
   expect_identical(getConditionallyIndependentSets(m, "x[1]", givenNodes = c("w[1]", "y"), unknownAsGiven=FALSE),
                    list(c("x[1]", "x[2]", "x[3]", "x[4]")))
   expect_identical(getConditionallyIndependentSets(m, givenNodes = c("y"), unknownAsGiven=FALSE),
-                   list(paste0("x[", 1:4, "]")))
+                   list(paste0("x[", c(2:4,1), "]")))
   expect_identical(getConditionallyIndependentSets(m, 'z[1]'), list())
   expect_true(nimble:::testConditionallyIndependentSets(m, getConditionallyIndependentSets(m)))
 
@@ -257,9 +257,9 @@ test_that("getConditionallyIndependentSets works in double-state state-space mod
   # expect_identical(getConditionallyIndependentSets(m, omit = "w[2]"),
   #                  list(c("x[2]", "x[3]", "w[3]", "x[4]", "w[4]")))
   expect_identical(getConditionallyIndependentSets(m, givenNodes = c("y", "w[3]"), unknownAsGiven=FALSE),
-                   list(c("x[1]", "w[1]", "x[2]", "x[3]", "x[4]", "w[2]", "w[4]")))
+                   list(c("x[2]", "x[3]", "x[4]", "x[1]", "w[2]", "w[4]", "w[1]")))
   expect_identical(getConditionallyIndependentSets(m, givenNodes = c("y", "x[3]", "w[3]"), unknownAsGiven=FALSE),
-                   list(c("x[1]", "w[1]", "x[2]", "w[2]"), c("x[4]", "w[4]")))
+                   list(c("x[2]", "x[1]", "w[2]", "w[1]"), c("x[4]", "w[4]")))
   expect_true(nimble:::testConditionallyIndependentSets(m, getConditionallyIndependentSets(m)))
 
   SMN <- setupMargNodes(m)
@@ -282,16 +282,16 @@ test_that("getConditionallyIndependentSets works in model with LHSinferred (aka
 
   expect_identical(getConditionallyIndependentSets(m3), list()) # Note there are no latent nodes
   expect_identical(getConditionallyIndependentSets(m3, "a[1:3]", givenNodes = "sig2", unknownAsGiven=FALSE),
-                   list(c("a[1:3]", "var", "c[1]", "c[2]")))
+                   list(c("a[1:3]", "c[2]", "c[1]", "var")))
   expect_identical(getConditionallyIndependentSets(m3, "a[1]", givenNodes = "sig2", unknownAsGiven=FALSE),
-                   list(c("a[1:3]", "var", "c[1]", "c[2]")))
+                   list(c("a[1:3]", "c[2]", "c[1]", "var")))
   expect_identical(getConditionallyIndependentSets(m3, "a[2]", givenNodes = "sig2", unknownAsGiven=FALSE),
-                   list(c("a[1:3]", "var", "c[1]", "c[2]")))
+                   list(c("a[1:3]", "c[2]", "c[1]", "var")))
 # Revisit
   #  expect_identical(getConditionallyIndependentSets(m3, "b[1]", givenNodes = "sig2", unknownAsGiven=FALSE),
 #                   list(c("var", "a[1:3]", "c[2]", "c[1]")))
   expect_identical(getConditionallyIndependentSets(m3, m3$getNodes(stochOnly = TRUE), givenNodes = "sig2"),
-                   list(c("a[1:3]" , "var", "c[1]", "c[2]"), c("a2")))
+                   list(c("a[1:3]" , "c[2]", "c[1]", "var"), c("a2")))
 })
 
 test_that("getConditionallyIndependentSets works for tweaked pump model", {
@@ -349,7 +349,7 @@ test_that("getConditionallyIndependentSets works with unknownAsGiven=TRUE or FAL
 
   expect_identical(
     getConditionallyIndependentSets(m, "a", givenNodes = c("y"), unknownAsGiven=FALSE),
-    list(c("mu", paste0("a[", 1:4, "]"))))
+    list(c(paste0("a[", 1:4, "]"), "mu")))
 
   expect_identical(
     getConditionallyIndependentSets(m, "a", givenNodes = c("y"), unknownAsGiven=TRUE),
@@ -491,7 +491,7 @@ test_that("`setupMargNodes` handling of missing/extra latents", {
     expect_message(result <- setupMargNodes(m, randomEffectsNodes = 'b[1]', calcNodes = c('b[1]','y')),
                    "they should be included")
     expect_identical(result$randomEffectsNodes[[1]]$toNodeChars(), c("b[1]"))
-    expect_identical(unlist(lapply(result$paramNodes, \(x) x$toNodeChars())), c("b[2]", "mu"))
+    expect_identical(unlist(lapply(result$paramNodes, \(x) x$toNodeChars())), c("mu", "b[2]"))
  
     expect_message(result <- setupMargNodes(m, randomEffectsNodes = 'b'), "they are not needed for the provided")
     expect_identical(result$randomEffectsNodes[[1]]$toNodeChars(), c("b[1]", "b[2]"))
@@ -524,12 +524,12 @@ test_that("getConditionallyIndependentSets works with grouped data", {
     
     expect_identical(getConditionallyIndependentSets(m, givenNodes = m$getNodes(c('mu0','y')),
                                                      nodes=m$getNodes(c('sigma','mu'))),
-                                        list(c("mu[1]","mu[2]","sigma"))) # {sigma,mu[1:2]}
+                                        list(c("sigma","mu[1]","mu[2]"))) # {sigma,mu[1:2]}
     
     expect_identical(getConditionallyIndependentSets(m, givenNodes = m$getNodes(c('mu0','sigma')),
                                                      nodes=m$getNodes(c('y','mu'))),
-                     list(c("mu[1]",paste0("y[",1:5,"]")),
-                          c("mu[2]", paste0("y[",6:10,"]"))))            # {mu[1],y[1],y[2:5]}, {mu[2],y[6],y[7:10]}
+                     list(c(paste0("y[",1:5,"]"), "mu[1]"),
+                          c(paste0("y[",6:10,"]"), "mu[2]")))            # {mu[1],y[1],y[2:5]}, {mu[2],y[6],y[7:10]}
     
     expect_identical(getConditionallyIndependentSets(m, givenNodes = m$getNodes(c('mu0','sigma')),
                                                      nodes=m$getNodes(c('mu','y'))),
@@ -571,8 +571,8 @@ test_that("getConditionallyIndependentSets works with grouped data", {
     result <- getConditionallyIndependentSets(m, givenNodes = m$getNodes(c('mu0','sigma')),
                                               nodes=m$getNodes(c('y','mu')))
     expect_identical(result,
-                                 list(c("mu[1]", paste0("y[",1:5,"]")),
-                                      c("mu[2]", paste0("y[",6:10,"]"))))  
+                                 list(c(paste0("y[",1:5,"]"), "mu[1]"),
+                                      c(paste0("y[",6:10,"]"), "mu[2]")))  
                                         # {mu[1],y[1],y[2:5]}, {mu[2],y[6],y[7:10]}
  
     # Now with intermediate deterministic nodes.
@@ -594,12 +594,12 @@ test_that("getConditionallyIndependentSets works with grouped data", {
     
     expect_identical(getConditionallyIndependentSets(m, givenNodes = m$getNodes(c('mu0','y')),
                                                      nodes=m$getNodes(c('sigma','mu'))),
-                                        list(c("mu[1]","mu[2]","sigma"))) # {sigma,mu[1:2]}
+                                        list(c("sigma", "mu[1]","mu[2]"))) # {sigma,mu[1:2]}
     
     expect_identical(getConditionallyIndependentSets(m, givenNodes = m$getNodes(c('mu0','sigma')),
                                                      nodes=m$getNodes(c('y','mu'))),
-                     list(c("mu[1]", paste0("y[",1:5,"]")),
-                          c("mu[2]", paste0("y[",6:10,"]"))))            # {mu[1],y[1],y[2:5]}, {mu[2],y[6],y[7:10]}
+                     list(c(paste0("y[",1:5,"]"), "mu[1]"),
+                          c(paste0("y[",6:10,"]"), "mu[2]")))            # {mu[1],y[1],y[2:5]}, {mu[2],y[6],y[7:10]}
     
     expect_identical(getConditionallyIndependentSets(m, givenNodes = m$getNodes(c('mu0','sigma')),
                                                      nodes=m$getNodes(c('mu','y'))),
diff --git a/nimbleModel/tests/testthat/test-varRange.R b/nimbleModel/tests/testthat/test-varRange.R
index 841a527..0d599e6 100644
--- a/nimbleModel/tests/testthat/test-varRange.R
+++ b/nimbleModel/tests/testthat/test-varRange.R
@@ -203,14 +203,14 @@ test_that("toVarChars works correctly", {
                 y[idx[k], i+1,3,j,i,2:4]~ dmnorm(z[1:3],pr[1:3,1:3])
     })
     
-    md <- modelDefClass$new(code, constants = list(idx = c(2,5,4)))
-    vr <- getDependencies(md, 'z', self=FALSE)[[1]]
-    expect_identical(vr$toVarChars(),
+    m <- nimbleModel(code, constants = list(idx = c(2,5,4)))
+    nr <- getDependencies(m, 'z', self=FALSE)[[1]]
+    expect_identical(nr$toVarChars(),
                      c("y[2, 2, 3, 1:2, 1, 2:4]", "y[4, 2, 3, 1:2, 1, 2:4]", "y[5, 2, 3, 1:2, 1, 2:4]", "y[2, 3, 3, 1:2, 2, 2:4]", "y[4, 3, 3, 1:2, 2, 2:4]", "y[5, 3, 3, 1:2, 2, 2:4]", "y[2, 4, 3, 1:2, 3, 2:4]", "y[4, 4, 3, 1:2, 3, 2:4]", "y[5, 4, 3, 1:2, 3, 2:4]", "y[2, 5, 3, 1:2, 4, 2:4]", "y[4, 5, 3, 1:2, 4, 2:4]", "y[5, 5, 3, 1:2, 4, 2:4]"))
 
-    tmp <- vr$toVarChars()
+    tmp <- nr$toVarChars()
     result <- unlist(lapply(tmp, function(x) varRangeClass$new(x)$toVarChars(expandScalars = TRUE)))
-    expect_identical(sort(vr$toVarChars(expandScalars = TRUE)), sort(result))
+    expect_identical(sort(nr$toVarChars(expandScalars = TRUE)), sort(result))
                      
 
     code <- quote({
@@ -218,16 +218,16 @@ test_that("toVarChars works correctly", {
             for(j in 1:3)
                 y[j,i] ~ dnorm(0,1)
     })
-    md <- modelDefClass$new(code)
-    vr <- getDependencies(md, 'y')[[1]]
+    m <- nimbleModel(code)
+    nr <- getDependencies(m, 'y')[[1]]
     gr <- expand.grid(1:3, 1:2)
-    expect_identical(vr$toVarChars(expandScalars = TRUE),
+    expect_identical(nr$toVarChars(expandScalars = TRUE),
                      paste0("y[", gr[,1], ", ", gr[,2], "]"))
 
     vr <- varRangeClass$new(list(newIndexRange(quote(1:2)),newIndexRange(quote(1:3))),
                             rangeToIndexSlot = c(2,1), varName = 'y')
     gr <- expand.grid(1:3,1:2)
-    expect_identical(vr$toVarChars(expandScalars = TRUE),
+    expect_identical(nr$toVarChars(expandScalars = TRUE),
                      paste0("y[", gr[,1], ", ", gr[,2], "]"))
     
 })