diff --git a/R/Information.R b/R/Information.R
index 41492d9be..39ea5f256 100644
--- a/R/Information.R
+++ b/R/Information.R
@@ -140,7 +140,7 @@ MeilaMutualInformation <- function(split1, split2) {
     jointEntropies
   
   # Return:
-  if (abs(mutualInformation) < .Machine$double.eps^0.5) 0 else mutualInformation
+  if (abs(mutualInformation) < .Machine[["double.eps"]]^0.5) 0 else mutualInformation
 }
 
 #' Variation of information for all split pairings
diff --git a/R/MSTSegments.R b/R/MSTSegments.R
index 22f584eeb..0610599d4 100644
--- a/R/MSTSegments.R
+++ b/R/MSTSegments.R
@@ -69,7 +69,7 @@ StrainCol <- function(distances, mapping, mstEnds = MSTEdges(distances),
   })
   strain <- logStrain - median(logStrain[is.finite(logStrain)])
   # Infinite values arise when orig == 0
-  maxVal <- max(abs(strain[is.finite(strain)])) + sqrt(.Machine$double.eps)
+  maxVal <- max(abs(strain[is.finite(strain)])) + sqrt(.Machine[["double.eps"]])
   nCols <- length(palette)
   bins <- cut(strain, seq(-maxVal, maxVal, length.out = nCols))
   
diff --git a/R/VisualizeMatching.R b/R/VisualizeMatching.R
index e242983ae..341dd25c4 100644
--- a/R/VisualizeMatching.R
+++ b/R/VisualizeMatching.R
@@ -46,14 +46,14 @@ VisualizeMatching <- function(Func, tree1, tree2, setPar = TRUE,
                               ...) {
   
   splits1 <- as.Splits(tree1)
-  edge1 <- tree1$edge
+  edge1 <- tree1[["edge"]]
   child1 <- edge1[, 2]
   nTip <- attr(splits1, "nTip")
   splitEdges1 <- vapply(as.integer(rownames(splits1)),
                         function(node) which(child1 == node), integer(1))
   
   splits2 <- as.Splits(tree2, tipLabels = tree1)
-  edge2 <- tree2$edge
+  edge2 <- tree2[["edge"]]
   child2 <- edge2[, 2]
   splitEdges2 <- vapply(as.integer(rownames(splits2)),
                         function(node) which(child2 == node), integer(1))
diff --git a/R/kmeanspp.R b/R/kmeanspp.R
index 0628de62b..e971d69d5 100644
--- a/R/kmeanspp.R
+++ b/R/kmeanspp.R
@@ -58,7 +58,7 @@ KMeansPP.matrix <- function(x, k = 2, nstart = 10, ...) {
     }
     
     proposal <- kmeans(x, centers = x[centres, ], ...)
-    if (proposal$tot.withinss < ret$tot.withinss){
+    if (proposal[["tot.withinss"]] < ret[["tot.withinss"]]){
       ret <- proposal
     }
   }
@@ -94,7 +94,7 @@ KMeansPP.dist <- function(x, k = 2, nstart = 10, ...) {
     }
     
     proposal <- kmeans(x, centers = d[centres, ], ...)
-    if (proposal$tot.withinss < ret$tot.withinss){
+    if (proposal[["tot.withinss"]] < ret[["tot.withinss"]]){
       ret <- proposal
     }
   }
diff --git a/R/plot.R b/R/plot.R
index 6d6149682..397039248 100644
--- a/R/plot.R
+++ b/R/plot.R
@@ -21,25 +21,25 @@ TreeDistPlot <- function(tr, title = NULL, bold = NULL, leaveRoom = FALSE,
                          prune = integer(0), graft = integer(0),
                          edge.color = "black", edge.width = NULL, ...) {
   
-  nEdge <- dim(tr$edge)[1]
-  if (is.null(tr$edge.length)) {
-    tr$edge.length <- rep(1, nEdge)
+  nEdge <- dim(tr[["edge"]])[1]
+  if (is.null(tr[["edge.length"]])) {
+    tr[["edge.length"]] <- rep(1, nEdge)
   }
   if (is.null(edge.width)) {
-    edge.width <- if (is.null(tr$edge.width)) {
+    edge.width <- if (is.null(tr[["edge.width"]])) {
       rep(1, nEdge)
     } else {
-      tr$edge.width
+      tr[["edge.width"]]
     }
   }
   if (length(edge.color) == 1) {
     edge.color <- rep(edge.color, nEdge)
   }
-  nTip <- length(tr$tip.label)
-  if (all(tr$tip.label %in% LETTERS)) {
-    tr$tip.label <- match(tr$tip.label, LETTERS)
-  } else if (all(tr$tip.label %in% letters)) {
-    tr$tip.label <- match(tr$tip.label, letters)
+  nTip <- length(tr[["tip.label"]])
+  if (all(tr[["tip.label"]] %in% LETTERS)) {
+    tr[["tip.label"]] <- match(tr[["tip.label"]], LETTERS)
+  } else if (all(tr[["tip.label"]] %in% letters)) {
+    tr[["tip.label"]] <- match(tr[["tip.label"]], letters)
   }
 
   if (length(prune) > 0 || length(graft) > 0) {
@@ -51,7 +51,7 @@ TreeDistPlot <- function(tr, title = NULL, bold = NULL, leaveRoom = FALSE,
                "#6DB6FF", "#B6DBFF", "#920000", "#924900", "#DB6D00", "#24FF24",
                "#FFFF6D") # Ternary::cbPalette15[-c(4, 7)]
 
-  tipNumbers <- tr$tip.label
+  tipNumbers <- tr[["tip.label"]]
   font <- rep(1L, length(tipNumbers))
   if (!is.null(bold)) {
     font[tipNumbers %in% bold] <- 4L
@@ -67,7 +67,7 @@ TreeDistPlot <- function(tr, title = NULL, bold = NULL, leaveRoom = FALSE,
                         warning("Leaves of `tr` must be labelled with integers")
                       })
   
-  tr$tip.label <- LETTERS[tipInts]
+  tr[["tip.label"]] <- LETTERS[tipInts]
   plot.phylo(tr, tip.color = tipCols[tipInts],
              main = title, cex.main = 0.8, font = font, 
              edge.width = edge.width, edge.color = edge.color, 
diff --git a/R/spectral_clustering.R b/R/spectral_clustering.R
index b08970311..b8faa6957 100644
--- a/R/spectral_clustering.R
+++ b/R/spectral_clustering.R
@@ -61,7 +61,7 @@ SpectralEigens <- function(D, nn = 10L, nEig = 2L) {
   ei <- eigen(L, symmetric = TRUE) # 3. Compute the eigenvectors and values of L
   
   # Return the eigenvectors of the n_eig smallest eigenvalues:
-  ei$vectors[, nrow(L) - rev(seq_len(nEig))]
+  ei[["vectors"]][, nrow(L) - rev(seq_len(nEig))]
 }
 
 #' @export
diff --git a/R/tree_distance.R b/R/tree_distance.R
index a3bd34277..bf8ded308 100644
--- a/R/tree_distance.R
+++ b/R/tree_distance.R
@@ -35,10 +35,10 @@ GeneralizedRF <- function(splits1, splits2, nTip, PairScorer,
   nSplits2 <- dim(splits2)[1]
   
   solution <- PairScorer(splits1, splits2, nTip, ...)
-  ret <- solution$score
+  ret <- solution[["score"]]
   
   if (reportMatching) {
-    matching <- solution$matching
+    matching <- solution[["matching"]]
     matching[matching > nSplits2 | matching == 0L] <- NA
     if (nSplits1 < nSplits2) {
       matching <- matching[seq_len(nSplits1)]
@@ -51,7 +51,7 @@ GeneralizedRF <- function(splits1, splits2, nTip, PairScorer,
       for (j in seq_len(nSplits2)) {
         pairScores[i, j] <- PairScorer(splits1[i, , drop = FALSE], 
                                        splits2[j, , drop = FALSE],
-                                       nTip = nTip, ...)$score
+                                       nTip = nTip, ...)[["score"]]
       }
     }
     attr(ret, "pairScores") <- pairScores
diff --git a/R/tree_distance_info.R b/R/tree_distance_info.R
index 183515dd7..85304df5d 100644
--- a/R/tree_distance_info.R
+++ b/R/tree_distance_info.R
@@ -221,7 +221,7 @@ DifferentPhylogeneticInfo <- function(tree1, tree2 = NULL, normalize = FALSE,
                        infoInBoth = treesIndependentInfo,
                        InfoInTree = SplitwiseInfo, Combine = "+")
   
-  ret[ret < .Machine$double.eps ^ 0.5] <- 0 # Catch floating point inaccuracy
+  ret[ret < .Machine[["double.eps"]] ^ 0.5] <- 0 # Catch floating point inaccuracy
   attributes(ret) <- attributes(spi)
   
   # Return:
@@ -246,7 +246,7 @@ ClusteringInfoDistance <- function(tree1, tree2 = NULL, normalize = FALSE,
                        infoInBoth = treesIndependentInfo,
                        InfoInTree = ClusteringEntropy, Combine = "+")
   
-  ret[ret < .Machine$double.eps ^ 0.5] <- 0 # Handle floating point inaccuracy
+  ret[ret < .Machine[["double.eps"]] ^ 0.5] <- 0 # Handle floating point inaccuracy
   attributes(ret) <- attributes(mci)
   
   # Return:
diff --git a/R/tree_distance_kendall-colijn.R b/R/tree_distance_kendall-colijn.R
index 3657c0828..283bd2ad4 100644
--- a/R/tree_distance_kendall-colijn.R
+++ b/R/tree_distance_kendall-colijn.R
@@ -82,8 +82,8 @@ KendallColijn <- function(tree1, tree2 = NULL, Vector = KCVector) {
   
   if (inherits(tree1, "phylo")) {
     if (inherits(tree2, "phylo")) {
-      if (length(tree1$tip.label) != length(tree2$tip.label) || 
-          length(setdiff(tree1$tip.label, tree2$tip.label)) > 0) {
+      if (length(tree1[["tip.label"]]) != length(tree2[["tip.label"]]) || 
+          length(setdiff(tree1[["tip.label"]], tree2[["tip.label"]])) > 0) {
         stop("Leaves must bear identical labels.")
       }
       .EuclideanDistance(Vector(tree1) - Vector(tree2))
@@ -92,18 +92,18 @@ KendallColijn <- function(tree1, tree2 = NULL, Vector = KCVector) {
         0
       } else {
         apply(Vector(tree1) - vapply(tree2, Vector,
-                                     FunValue(length(tree1$tip.label))),
+                                     FunValue(length(tree1[["tip.label"]]))),
               2L, .EuclideanDistance)
       }
     }
   } else {
     if (inherits(tree2, "phylo")) {
       apply(Vector(tree2) - vapply(tree1, Vector,
-                                   FunValue(length(tree2$tip.label))),
+                                   FunValue(length(tree2[["tip.label"]]))),
             2L, .EuclideanDistance)
     } else if (is.null(tree2)) {
       
-      treeVec <- vapply(tree1, Vector, FunValue(length(tree1[[1]]$tip.label)))
+      treeVec <- vapply(tree1, Vector, FunValue(length(tree1[[1]][["tip.label"]])))
       nTree <- length(tree1)
       ret <- matrix(0, nTree, nTree)
       is <- combn(seq_len(nTree), 2)
@@ -116,8 +116,8 @@ KendallColijn <- function(tree1, tree2 = NULL, Vector = KCVector) {
       ret
     }
       else {
-      vector1 <- vapply(tree1, Vector, FunValue(length(tree1[[1]]$tip.label)))
-      vector2 <- vapply(tree2, Vector, FunValue(length(tree2[[1]]$tip.label)))
+      vector1 <- vapply(tree1, Vector, FunValue(length(tree1[[1]][["tip.label"]])))
+      vector2 <- vapply(tree2, Vector, FunValue(length(tree2[[1]][["tip.label"]])))
       apply(vector2, 2, function(i) 
         apply(vector1, 2, function(j) 
           .EuclideanDistance(i - j)))
@@ -135,12 +135,12 @@ KendallColijn <- function(tree1, tree2 = NULL, Vector = KCVector) {
 #' @export
 KCVector <- function(tree) {
   tree <- Preorder(tree)
-  edge <- tree$edge
+  edge <- tree[["edge"]]
   parent <- edge[, 1L]
   child <- edge[, 2L]
   root <- parent[1]
   nTip <- root - 1L
-  tipOrder <- order(tree$tip.label)
+  tipOrder <- order(tree[["tip.label"]])
   is <- combn(tipOrder, 2)
   
   ancestors <- AllAncestors(parent, child)
@@ -163,13 +163,13 @@ PathVector <- function(tree) {
     stop("`tree` must be of class `phylo`")
   }
   tree <- Preorder(tree)
-  edge <- tree$edge
+  edge <- tree[["edge"]]
   parent <- edge[, 1L]
   child <- edge[, 2L]
   root <- parent[1]
   nTip <- root - 1L
   tipAncs <- seq_len(nTip)
-  tipOrder <- order(tree$tip.label)
+  tipOrder <- order(tree[["tip.label"]])
   is <- combn(tipOrder, 2)
   
   ancestors <- AllAncestors(parent, child)
@@ -192,7 +192,7 @@ PathVector <- function(tree) {
 #' @importFrom TreeTools as.Splits
 #' @export
 SplitVector <- function(tree) {
-  tipLabel <- tree$tip.label
+  tipLabel <- tree[["tip.label"]]
   nTip <- length(tipLabel)
   splits <- as.logical(as.Splits(tree, tipLabel[order(tipLabel)]))
   splits <- rbind(splits, !splits)
@@ -245,4 +245,4 @@ KCDiameter.list <- function(tree) {
 #' @export
 KCDiameter.multiPhylo <- function(tree) {
   vapply(tree, KCDiameter, double(1))
-}
\ No newline at end of file
+}
diff --git a/R/tree_distance_mast.R b/R/tree_distance_mast.R
index 77ee4653e..4cdd06e01 100644
--- a/R/tree_distance_mast.R
+++ b/R/tree_distance_mast.R
@@ -55,32 +55,32 @@ MASTSize <- function(tree1, tree2 = tree1, rooted = TRUE) {
 #' @importFrom ape drop.tip
 #' @importFrom TreeTools Postorder RenumberTips TreeIsRooted RootOnNode
 .MASTSizeSingle <- function(tree1, tree2, rooted = TRUE,
-                             tipLabels = tree1$tip.label,
+                             tipLabels = tree1[["tip.label"]],
                              ...) {
   label1 <- tipLabels
-  label2 <- tree2$tip.label
+  label2 <- tree2[["tip.label"]]
   
   tree1 <- drop.tip(tree1, setdiff(label1, label2))
-  label1 <- tree1$tip.label
+  label1 <- tree1[["tip.label"]]
   tree2 <- RenumberTips(drop.tip(tree2, setdiff(label2, label1)), label1)
   
   nTip <- length(label1)
   
   if (!rooted) {
     if (!TreeIsRooted(tree1)) {
-      tree1 <- RootOnNode(tree1, node = tree1$edge[nTip + nTip - 2L], TRUE)
+      tree1 <- RootOnNode(tree1, node = tree1[["edge"]][nTip + nTip - 2L], TRUE)
     }
-    postorderEdge1 <- Postorder(tree1$edge)
+    postorderEdge1 <- Postorder(tree1[["edge"]])
     tree2 <- Preorder(tree2)
-    max(vapply(tree2$edge[, 2], function(node)
+    max(vapply(tree2[["edge"]][, 2], function(node)
       .MASTSizeEdges(postorderEdge1,
-                     RootOnNode(tree2, node = node, TRUE)$edge,
+                     RootOnNode(tree2, node = node, TRUE)[["edge"]],
                      nTip = nTip), 0L))
   } else {
     if (!TreeIsRooted(tree1) || !TreeIsRooted(tree2)) {
       stop("Both trees must be rooted if rooted = TRUE")
     }
-    .MASTSizeEdges(Postorder(tree1$edge), tree2$edge, nTip)
+    .MASTSizeEdges(Postorder(tree1[["edge"]]), tree2[["edge"]], nTip)
   }
 }
 
diff --git a/R/tree_distance_msi.R b/R/tree_distance_msi.R
index 593430c71..eab4bdacc 100644
--- a/R/tree_distance_msi.R
+++ b/R/tree_distance_msi.R
@@ -29,7 +29,7 @@ MatchingSplitInfoDistance <- function(tree1, tree2 = NULL,
                        infoInBoth = treesIndependentInfo,
                        InfoInTree = SplitwiseInfo, Combine = "+")
   
-  ret[ret < .Machine$double.eps^0.5] <- 0 # In case of floating point inaccuracy
+  ret[ret < .Machine[["double.eps"]]^0.5] <- 0 # In case of floating point inaccuracy
   attributes(ret) <- attributes(msi)
   # Return:
   ret
diff --git a/R/tree_distance_nni.R b/R/tree_distance_nni.R
index b96857d36..b9c69b96a 100644
--- a/R/tree_distance_nni.R
+++ b/R/tree_distance_nni.R
@@ -73,10 +73,10 @@ NNIDist <- function(tree1, tree2 = tree1) {
 #' @importFrom TreeTools Postorder RenumberTips
 #' @importFrom ape Nnode.phylo
 .NNIDistSingle <- function(tree1, tree2, nTip, ...) {
-  tree2 <- RenumberTips(tree2, tree1$tip.label)
+  tree2 <- RenumberTips(tree2, tree1[["tip.label"]])
   
-  edge1 <- Postorder(tree1$edge)
-  edge2 <- Postorder(tree2$edge)
+  edge1 <- Postorder(tree1[["edge"]])
+  edge2 <- Postorder(tree2[["edge"]])
   
   cpp_nni_distance(edge1, edge2, nTip)
 }
diff --git a/R/tree_distance_rf.R b/R/tree_distance_rf.R
index 6205e648d..02721e6ef 100644
--- a/R/tree_distance_rf.R
+++ b/R/tree_distance_rf.R
@@ -74,7 +74,7 @@ InfoRobinsonFoulds <- function(tree1, tree2 = NULL, similarity = FALSE,
   }
   
   # In case of floating point inaccuracy
-  unnormalized[unnormalized < .Machine$double.eps^0.5] <- 0
+  unnormalized[unnormalized < .Machine[["double.eps"]]^0.5] <- 0
   
   # Return:
   NormalizeInfo(unnormalized, tree1, tree2, how = normalize,
diff --git a/R/tree_information.R b/R/tree_information.R
index 0c747268f..08fab0a95 100644
--- a/R/tree_information.R
+++ b/R/tree_information.R
@@ -136,7 +136,7 @@ SplitwiseInfo <- function(x, p = NULL, sum = TRUE) {
     if (p == FALSE) {
       p <- NULL
     } else {
-      np <- tree$node.label[as.integer(names(splits)) - NTip(tree)]
+      np <- tree[["node.label"]][as.integer(names(splits)) - NTip(tree)]
       if (is.null(np)) {
         np <- rep_len(p, length(splits))
       }