diff --git a/bct/algorithms/clustering.py b/bct/algorithms/clustering.py
index 484cc23..5732b3d 100644
--- a/bct/algorithms/clustering.py
+++ b/bct/algorithms/clustering.py
@@ -250,7 +250,7 @@ def clustering_coef_wu_sign(W, coef_type='default'):
     '''
     Returns the weighted clustering coefficient generalized or separated
     for positive and negative weights.
-  
+
     Three Algorithms are supported; herefore referred to as default, zhang,
     and costantini.
 
@@ -384,9 +384,8 @@ def consensus_und(D, tau, reps=1000, seed=None):
     reps : int
         number of times the clustering algorithm is reapplied. default value
         is 1000.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -485,7 +484,7 @@ def get_components(A, no_depend=False):
     if not np.all(A == A.T):  # ensure matrix is undirected
         raise BCTParamError('get_components can only be computed for undirected'
                             ' matrices.  If your matrix is noisy, correct it with np.around')
-    
+
     A = binarize(A, copy=True)
     n = len(A)
     np.fill_diagonal(A, 1)
@@ -503,7 +502,7 @@ def get_components(A, no_depend=False):
         temp.append(item)
         union_sets = temp
 
-    comps = np.array([i+1 for v in range(n) for i in 
+    comps = np.array([i+1 for v in range(n) for i in
         range(len(union_sets)) if v in union_sets[i]])
     comp_sizes = np.array([len(s) for s in union_sets])
 
diff --git a/bct/algorithms/core.py b/bct/algorithms/core.py
index 423f334..fb0dcef 100644
--- a/bct/algorithms/core.py
+++ b/bct/algorithms/core.py
@@ -166,9 +166,8 @@ def core_periphery_dir(W, gamma=1, C0=None, seed=None):
         0 < gamma < 1 detects large core, small periphery
     C0 : NxN np.ndarray
         Initial core structure
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
     '''
     rng = get_rng(seed)
     n = len(W)
diff --git a/bct/algorithms/generative.py b/bct/algorithms/generative.py
index fe5c805..f5b46ed 100644
--- a/bct/algorithms/generative.py
+++ b/bct/algorithms/generative.py
@@ -10,7 +10,7 @@
 
 
 @due.dcite(BibTeX(BETZEL2016), description="Generative models")
-def generative_model(A, D, m, eta, gamma=None, model_type='matching', 
+def generative_model(A, D, m, eta, gamma=None, model_type='matching',
     model_var='powerlaw', epsilon=1e-6, copy=True, seed=None):
     '''
     Generates synthetic networks using the models described in
@@ -32,10 +32,10 @@ def generative_model(A, D, m, eta, gamma=None, model_type='matching',
     D : np.ndarray
         Matrix of euclidean distances or other distances between nodes
     m : int
-        Number of connections that should be present in the final synthetic 
+        Number of connections that should be present in the final synthetic
         network
     eta : np.ndarray
-        A vector describing a range of values to estimate for eta, the 
+        A vector describing a range of values to estimate for eta, the
         hyperparameter describing exponential weighting of the euclidean
         distance.
     gamma : np.ndarray
@@ -44,7 +44,7 @@ def generative_model(A, D, m, eta, gamma=None, model_type='matching',
         algorithm. If model_type='euclidean' or another distance metric,
         this can be None.
     model_type : Enum(str)
-        euclidean : Uses only euclidean distances to generate connection 
+        euclidean : Uses only euclidean distances to generate connection
             probabilities
         neighbors : count of common neighbors
         matching : matching index, the normalized overlap in neighborhoods
@@ -67,16 +67,15 @@ def generative_model(A, D, m, eta, gamma=None, model_type='matching',
     copy : bool
         Some algorithms add edges directly to the input matrix. Set this flag
         to make a copy of the input matrix instead. Defaults to True.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
     '''
     rng = get_rng(seed)
     if copy:
         A = A.copy()
 
     n = len(D)
-    
+
     #These parameters don't do any of the voronoi narrowing.
     #Its a list of eta values paired with gamma values.
     #To try 3 eta and 3 gamma pairs, should use 9 list values.
@@ -93,7 +92,7 @@ def k_avg(K):
             epsilon)
 
     def k_diff(K):
-        return np.abs(np.tile(K, (n, 1)) - 
+        return np.abs(np.tile(K, (n, 1)) -
                       np.transpose(np.tile(K, (n, 1)))) + epsilon
 
     def k_max(K):
@@ -117,7 +116,7 @@ def s_diff(K, sc):
 
     def s_min(K, sc):
         return np.where(K < sc, K + epsilon, sc + epsilon)
-    
+
     def s_max(K, sc):
         #return np.max((K, sc.T), axis=0)
         return np.where(K > sc, K + epsilon, sc + epsilon)
@@ -169,12 +168,12 @@ def clu_gen(A, K, D, m, eta, gamma, model_var, x_fun):
         if mv1 in ('powerlaw', 'power_law'):
             Fd = D**eta
         elif mv1 in ('exponential',):
-            Fd = np.exp(eta*D) 
+            Fd = np.exp(eta*D)
 
         if mv2 in ('powerlaw', 'power_law'):
             Fk = K**gamma
         elif mv2 in ('exponential',):
-            Fk = np.exp(gamma*K) 
+            Fk = np.exp(gamma*K)
 
         c = clustering_coef_bu(A)
         k = np.sum(A, axis=1)
@@ -204,7 +203,7 @@ def clu_gen(A, K, D, m, eta, gamma, model_var, x_fun):
             c[k<=1] = 0
             bth[uu] = 1
             bth[vv] = 1
-    
+
             k_result = x_fun(c, bth)
 
             #print(np.shape(k_result))
@@ -239,12 +238,12 @@ def deg_gen(A, K, D, m, eta, gamma, model_var, s_fun):
         if mv1 in ('powerlaw', 'power_law'):
             Fd = D**eta
         elif mv1 in ('exponential',):
-            Fd = np.exp(eta*D) 
+            Fd = np.exp(eta*D)
 
         if mv2 in ('powerlaw', 'power_law'):
             Fk = K**gamma
         elif mv2 in ('exponential',):
-            Fk = np.exp(gamma*K) 
+            Fk = np.exp(gamma*K)
 
         P = Fd * Fk * np.logical_not(A)
         u,v = np.where(np.triu(np.ones((n,n)), 1))
@@ -258,7 +257,7 @@ def deg_gen(A, K, D, m, eta, gamma, model_var, s_fun):
 #        print(np.shape(np.where(A[u,v])), 'sqishy')
 #        print(np.shape(P), 'squnnaq')
 
-        #b[:mseed] = np.where(A[np.ix_(u,v)]) 
+        #b[:mseed] = np.where(A[np.ix_(u,v)])
         b[:mseed] = np.squeeze(np.where(A[u,v]))
         #print(mseed, m)
         for i in range(mseed, m):
@@ -318,16 +317,16 @@ def matching_gen(A, K, D, m, eta, gamma, model_var):
         if mv1 in ('powerlaw', 'power_law'):
             Fd = D**eta
         elif mv1 in ('exponential',):
-            Fd = np.exp(eta*D) 
+            Fd = np.exp(eta*D)
 
         if mv2 in ('powerlaw', 'power_law'):
             Fk = K**gamma
         elif mv2 in ('exponential',):
-            Fk = np.exp(gamma*K) 
+            Fk = np.exp(gamma*K)
 
         Ff = Fd * Fk * np.logical_not(A)
         u,v = np.where(np.triu(np.ones((n,n)), 1))
-    
+
         for ii in range(mseed, m):
             C = np.append(0, np.cumsum(Ff[u,v]))
             r = np.sum(rng.random_sample()*C[-1] >= C)
@@ -343,7 +342,7 @@ def matching_gen(A, K, D, m, eta, gamma, model_var):
             for i in range(len(updateuu)):
                 j = updateuu[i]
                 c2 = np.append(A[:,j], A[j,:])
-    
+
                 use = np.logical_or(c1, c2)
                 use[uu] = use[uu+n] = use[j] = use[j+n] = 0
                 ncon = np.sum(c1[use]) + np.sum(c2[use])
@@ -356,12 +355,12 @@ def matching_gen(A, K, D, m, eta, gamma, model_var):
             updatevv, = np.where(np.inner(A, A[:,vv]))
             np.delete(updatevv, np.where(updatevv == uu))
             np.delete(updatevv, np.where(updatevv == vv))
-        
+
             c1 = np.append(A[:,vv], A[vv,:])
             for i in range(len(updatevv)):
                 j = updatevv[i]
                 c2 = np.append(A[:,j], A[j,:])
-    
+
                 use = np.logical_or(c1, c2)
                 use[vv] = use[vv+n] = use[j] = use[j+n] = 0
                 ncon = np.sum(c1[use]) + np.sum(c2[use])
@@ -374,7 +373,7 @@ def matching_gen(A, K, D, m, eta, gamma, model_var):
             Ff = Fd * Fk * np.logical_not(A)
 
         return A
-    
+
     def neighbors_gen(A, K, D, m, eta, gamma, model_var):
         K += epsilon
 
@@ -388,16 +387,16 @@ def neighbors_gen(A, K, D, m, eta, gamma, model_var):
         if mv1 in ('powerlaw', 'power_law'):
             Fd = D**eta
         elif mv1 in ('exponential',):
-            Fd = np.exp(eta*D) 
+            Fd = np.exp(eta*D)
 
         if mv2 in ('powerlaw', 'power_law'):
             Fk = K**gamma
         elif mv2 in ('exponential',):
-            Fk = np.exp(gamma*K) 
+            Fk = np.exp(gamma*K)
 
         Ff = Fd * Fk * np.logical_not(A)
         u,v = np.where(np.triu(np.ones((n,n)), 1))
-    
+
         for ii in range(mseed, m):
             C = np.append(0, np.cumsum(Ff[u,v]))
             r = np.sum(rng.random_sample()*C[-1] >= C)
@@ -407,7 +406,7 @@ def neighbors_gen(A, K, D, m, eta, gamma, model_var):
 
             x = A[uu, :].astype(int)
             y = A[:, vv].astype(int)
-    
+
             K[uu, y] += 1
             K[y, uu] += 1
             K[vv, x] += 1
@@ -416,7 +415,7 @@ def neighbors_gen(A, K, D, m, eta, gamma, model_var):
             if mv2 in ('powerlaw', 'power_law'):
                 Fk = K**gamma
             elif mv2 in ('exponential',):
-                Fk = np.exp(gamma*K) 
+                Fk = np.exp(gamma*K)
 
             if mv2 in ('powerlaw', 'power_law'):
                 Ff[uu, y] = Ff[y, uu] = Fd[uu, y] * (K[uu, y] ** gamma)
@@ -474,12 +473,12 @@ def euclidean_gen(A, D, m, eta, model_var):
     elif model_type in ('clu-max', 'clu_max'):
         Kseed = k_max(clustering_coef_bu(A))
         for j, (ep, gp) in enumerate(zip(eta, gamma)):
-            B[:,:,j] = clu_gen(A, Kseed, D, m, ep, gp, model_var, x_max) 
+            B[:,:,j] = clu_gen(A, Kseed, D, m, ep, gp, model_var, x_max)
 
     elif model_type in ('clu-min', 'clu_min'):
         Kseed = k_min(clustering_coef_bu(A))
         for j, (ep, gp) in enumerate(zip(eta, gamma)):
-            B[:,:,j] = clu_gen(A, Kseed, D, m, ep, gp, model_var, x_min) 
+            B[:,:,j] = clu_gen(A, Kseed, D, m, ep, gp, model_var, x_min)
 
     elif model_type in ('clu-prod', 'clu_prod'):
         Kseed = k_prod(clustering_coef_bu(A))
@@ -495,7 +494,7 @@ def euclidean_gen(A, D, m, eta, model_var):
         Kseed = k_diff(np.sum(A, axis=1))
         for j, (ep, gp) in enumerate(zip(eta, gamma)):
             B[:,:,j] = deg_gen(A, Kseed, D, m, ep, gp, model_var, s_diff)
-    
+
     elif model_type in ('deg-max', 'deg_max'):
         Kseed = k_max(np.sum(A, axis=1))
         for j, (ep, gp) in enumerate(zip(eta, gamma)):
@@ -525,11 +524,11 @@ def euclidean_gen(A, D, m, eta, model_var):
 
     elif model_type in ('spatial', 'geometric', 'euclidean'):
         for j, ep in enumerate(eta):
-            B[:,:,j] = euclidean_gen(A, D, m, ep, model_var) 
+            B[:,:,j] = euclidean_gen(A, D, m, ep, model_var)
 
     return np.squeeze(B)
 
-def evaluate_generative_model(A, Atgt, D, eta, gamma=None, 
+def evaluate_generative_model(A, Atgt, D, eta, gamma=None,
     model_type='matching', model_var='powerlaw', epsilon=1e-6, seed=None):
     '''
     Generates synthetic networks with parameters provided and evaluates their
@@ -537,7 +536,7 @@ def evaluate_generative_model(A, Atgt, D, eta, gamma=None,
     Basically it takes the Kolmogorov-Smirnov statistics of 4 network
     measures; comparing the degree distributions, clustering coefficients,
     betweenness centrality, and Euclidean distances between connected regions.
-    
+
     The energy is globally low if the synthetic network matches the target.
     Energy is defined as the maximum difference across the four statistics.
     '''
@@ -548,11 +547,11 @@ def evaluate_generative_model(A, Atgt, D, eta, gamma=None,
     xb = betweenness_bin(Atgt)
     xe = D[np.triu(Atgt, 1) > 0]
 
-    B = generative_model(A, D, m, eta, gamma, model_type=model_type, 
+    B = generative_model(A, D, m, eta, gamma, model_type=model_type,
                          model_var=model_var, epsilon=epsilon, copy=True, seed=seed)
 
     #if eta != gamma then an error is thrown within generative model
-    
+
     nB = len(eta)
 
     if nB == 1:
@@ -562,7 +561,7 @@ def evaluate_generative_model(A, Atgt, D, eta, gamma=None,
 
     def kstats(x, y):
         bin_edges = np.concatenate([[-np.inf],
-                                    np.sort(np.concatenate((x, y))), 
+                                    np.sort(np.concatenate((x, y))),
                                     [np.inf]])
 
         bin_x,_ = np.histogram(x, bin_edges)
diff --git a/bct/algorithms/modularity.py b/bct/algorithms/modularity.py
index 6712377..efaec38 100644
--- a/bct/algorithms/modularity.py
+++ b/bct/algorithms/modularity.py
@@ -98,9 +98,8 @@ def community_louvain(W, gamma=1, ci=None, B='modularity', seed=None):
             'negative_sym' symmetric treatment of negative weights
             'negative_asym' asymmetric treatment of negative weights
         The default value is to use the Q-metric
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -607,9 +606,8 @@ def modularity_finetune_dir(W, ci=None, gamma=1, seed=None):
     gamma : float
         resolution parameter. default value=1. Values 0 <= gamma < 1 detect
         larger modules while gamma > 1 detects smaller modules.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -711,9 +709,8 @@ def modularity_finetune_und(W, ci=None, gamma=1, seed=None):
     gamma : float
         resolution parameter. default value=1. Values 0 <= gamma < 1 detect
         larger modules while gamma > 1 detects smaller modules.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -815,9 +812,8 @@ def modularity_finetune_und_sign(W, qtype='sta', gamma=1, ci=None, seed=None):
         larger modules while gamma > 1 detects smaller modules.
     ci : Nx1 np.ndarray | None
         initial community affiliation vector
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -949,9 +945,8 @@ def modularity_louvain_dir(W, gamma=1, hierarchy=False, seed=None):
         larger modules while gamma > 1 detects smaller modules.
     hierarchy : bool
         Enables hierarchical output. Defalut value=False
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1082,9 +1077,8 @@ def modularity_louvain_und(W, gamma=1, hierarchy=False, seed=None):
         larger modules while gamma > 1 detects smaller modules.
     hierarchy : bool
         Enables hierarchical output. Defalut value=False
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1222,9 +1216,8 @@ def modularity_louvain_und_sign(W, gamma=1, qtype='sta', seed=None):
     gamma : float
         resolution parameter. default value=1. Values 0 <= gamma < 1 detect
         larger modules while gamma > 1 detects smaller modules.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1392,9 +1385,8 @@ def modularity_probtune_und_sign(W, qtype='sta', gamma=1, ci=None, p=.45,
         initial community affiliation vector
     p : float
         probability of random node moves. Default value = 0.45
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
diff --git a/bct/algorithms/physical_connectivity.py b/bct/algorithms/physical_connectivity.py
index e97864b..aecaf27 100644
--- a/bct/algorithms/physical_connectivity.py
+++ b/bct/algorithms/physical_connectivity.py
@@ -93,9 +93,8 @@ def rentian_scaling(A, xyz, n, seed=None):
     n : int
         Number of partitions to compute. Each partition is a data point; you
         want a large enough number to adequately compute Rent's exponent.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
diff --git a/bct/algorithms/reference.py b/bct/algorithms/reference.py
index 6e33510..77876ac 100644
--- a/bct/algorithms/reference.py
+++ b/bct/algorithms/reference.py
@@ -27,9 +27,8 @@ def latmio_dir_connected(R, itr, D=None, seed=None):
     D : np.ndarray | None
         distance-to-diagonal matrix. Defaults to the actual distance matrix
         if not specified.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -149,9 +148,8 @@ def latmio_dir(R, itr, D=None, seed=None):
     D : np.ndarray | None
         distance-to-diagonal matrix. Defaults to the actual distance matrix
         if not specified.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -249,9 +247,8 @@ def latmio_und_connected(R, itr, D=None, seed=None):
     D : np.ndarray | None
         distance-to-diagonal matrix. Defaults to the actual distance matrix
         if not specified.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -384,9 +381,8 @@ def latmio_und(R, itr, D=None, seed=None):
     D : np.ndarray | None
         distance-to-diagonal matrix. Defaults to the actual distance matrix
         if not specified.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -488,9 +484,8 @@ def makeevenCIJ(n, k, sz_cl, seed=None):
         number of edges
     sz_cl : int
         size of clusters (must be power of 2)
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -566,9 +561,8 @@ def makefractalCIJ(mx_lvl, E, sz_cl, seed=None):
         connection density fall off per level
     sz_cl : int
         size of clusters (must be power of 2)
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -622,9 +616,8 @@ def makerandCIJdegreesfixed(inv, outv, seed=None):
         in-degree vector
     outv : Nx1 np.ndarray
         out-degree vector
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -706,9 +699,8 @@ def makerandCIJ_dir(n, k, seed=None):
         number of vertices
     K : int
         number of edges
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -738,9 +730,8 @@ def makerandCIJ_und(n, k, seed=None):
         number of vertices
     K : int
         number of edges
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -771,9 +762,8 @@ def makeringlatticeCIJ(n, k, seed=None):
         number of vertices
     K : int
         number of edges
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -829,9 +819,8 @@ def maketoeplitzCIJ(n, k, s, seed=None):
         number of edges
     s : float
         standard deviation of toeplitz
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -881,9 +870,8 @@ def null_model_dir_sign(W, bin_swaps=5, wei_freq=.1, seed=None):
         wei_freq == 0.1 implies that weights sorted each 10th step (faster,
             default value)
         wei_freq == 0 implies no sorting of weights (not recommended)
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1005,9 +993,8 @@ def null_model_und_sign(W, bin_swaps=5, wei_freq=.1, seed=None):
         wei_freq == 0.1 implies that weights sorted each 10th step (faster,
             default value)
         wei_freq == 0 implies no sorting of weights (not recommended)
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1131,9 +1118,8 @@ def randmio_dir_connected(R, itr, seed=None):
         directed binary/weighted connection matrix
     itr : int
         rewiring parameter. Each edge is rewired approximately itr times.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1223,9 +1209,8 @@ def randmio_dir(R, itr, seed=None):
         directed binary/weighted connection matrix
     itr : int
         rewiring parameter. Each edge is rewired approximately itr times.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1300,9 +1285,8 @@ def randmio_und_connected(R, itr, seed=None):
         undirected binary/weighted connection matrix
     itr : int
         rewiring parameter. Each edge is rewired approximately itr times.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1412,9 +1396,8 @@ def randmio_dir_signed(R, itr, seed=None):
         directed binary/weighted connection matrix
     itr : int
         rewiring parameter. Each edge is rewired approximately itr times.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1488,9 +1471,8 @@ def randmio_und(R, itr, seed=None):
         undirected binary/weighted connection matrix
     itr : int
         rewiring parameter. Each edge is rewired approximately itr times.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1571,9 +1553,8 @@ def randmio_und_signed(R, itr, seed=None):
         undirected binary/weighted connection matrix
     itr : int
         rewiring parameter. Each edge is rewired approximately itr times.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1636,9 +1617,8 @@ def randomize_graph_partial_und(A, B, maxswap, seed=None):
         mask; edges to avoid
     maxswap : int
         number of rewirings
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1711,9 +1691,8 @@ def randomizer_bin_und(R, alpha, seed=None):
         binary undirected connection matrix
     alpha : float
         fraction of edges to rewire
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
diff --git a/bct/nbs.py b/bct/nbs.py
index ce5006c..be397fd 100644
--- a/bct/nbs.py
+++ b/bct/nbs.py
@@ -27,7 +27,7 @@ def nbs_bct(x, y, thresh, k=1000, tail='both', paired=False, verbose=False, seed
     thresh : float
         minimum t-value used as threshold
     k : int
-        number of permutations used to estimate the empirical null 
+        number of permutations used to estimate the empirical null
         distribution
     tail : {'left', 'right', 'both'}
         enables specification of particular alternative hypothesis
@@ -39,9 +39,8 @@ def nbs_bct(x, y, thresh, k=1000, tail='both', paired=False, verbose=False, seed
         subject populations to have equal N. default value = False
     verbose : bool
         print some extra information each iteration. defaults value = False
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -55,42 +54,42 @@ def nbs_bct(x, y, thresh, k=1000, tail='both', paired=False, verbose=False, seed
         an adjacency matrix identifying the edges comprising each component.
         edges are assigned indexed values.
     null : Kx1 np.ndarray
-        A vector of K sampled from the null distribution of maximal component 
+        A vector of K sampled from the null distribution of maximal component
         size.
 
     Notes
     -----
-    ALGORITHM DESCRIPTION 
-    The NBS is a nonparametric statistical test used to isolate the 
-    components of an N x N undirected connectivity matrix that differ 
-    significantly between two distinct populations. Each element of the 
-    connectivity matrix stores a connectivity value and each member of 
-    the two populations possesses a distinct connectivity matrix. A 
-    component of a connectivity matrix is defined as a set of 
-    interconnected edges. 
-
-    The NBS is essentially a procedure to control the family-wise error 
-    rate, in the weak sense, when the null hypothesis is tested 
+    ALGORITHM DESCRIPTION
+    The NBS is a nonparametric statistical test used to isolate the
+    components of an N x N undirected connectivity matrix that differ
+    significantly between two distinct populations. Each element of the
+    connectivity matrix stores a connectivity value and each member of
+    the two populations possesses a distinct connectivity matrix. A
+    component of a connectivity matrix is defined as a set of
+    interconnected edges.
+
+    The NBS is essentially a procedure to control the family-wise error
+    rate, in the weak sense, when the null hypothesis is tested
     independently at each of the N(N-1)/2 edges comprising the undirected
-    connectivity matrix. The NBS can provide greater statistical power 
-    than conventional procedures for controlling the family-wise error 
+    connectivity matrix. The NBS can provide greater statistical power
+    than conventional procedures for controlling the family-wise error
     rate, such as the false discovery rate, if the set of edges at which
     the null hypothesis is rejected constitues a large component or
     components.
     The NBS comprises fours steps:
     1. Perform a two-sample T-test at each edge indepedently to test the
        hypothesis that the value of connectivity between the two
-       populations come from distributions with equal means. 
+       populations come from distributions with equal means.
     2. Threshold the T-statistic available at each edge to form a set of
-       suprathreshold edges. 
+       suprathreshold edges.
     3. Identify any components in the adjacency matrix defined by the set
-       of suprathreshold edges. These are referred to as observed 
-       components. Compute the size of each observed component 
-       identified; that is, the number of edges it comprises. 
+       of suprathreshold edges. These are referred to as observed
+       components. Compute the size of each observed component
+       identified; that is, the number of edges it comprises.
     4. Repeat K times steps 1-3, each time randomly permuting members of
-       the two populations and storing the size of the largest component 
+       the two populations and storing the size of the largest component
        identified for each permuation. This yields an empirical estimate
-       of the null distribution of maximal component size. A corrected 
+       of the null distribution of maximal component size. A corrected
        p-value for each observed component is then calculated using this
        null distribution.