diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 19fcdb8e..96fdd2b1 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -6,6 +6,7 @@ Version v3.1.0
 
 Improvements
 ~~~~~~~~~~~~
+- Support for ``numpy>=2``
 - Made ``Edges`` and ``EdgePopulation`` ``get`` functions more consistent
 
   - Both now return ``self.ids(query)`` if ``properties=None``
diff --git a/bluepysnap/circuit_validation.py b/bluepysnap/circuit_validation.py
index c40443f7..9f153a20 100644
--- a/bluepysnap/circuit_validation.py
+++ b/bluepysnap/circuit_validation.py
@@ -4,6 +4,7 @@
 """
 
 import logging
+from contextlib import contextmanager
 from pathlib import Path
 
 import h5py
@@ -21,6 +22,13 @@
 MAX_MISSING_FILES_DISPLAY = 10
 
 
+@contextmanager
+def _np_printoptions():
+    """Get numpy print options."""
+    with np.printoptions(legacy=("1.25" if np.version.version >= "2.0.0" else None)):
+        yield
+
+
 def _check_partial_circuit_config(config):
     return config.get("metadata", {}).get("status") == "partial"
 
@@ -347,11 +355,12 @@ def _check_edges_node_ids(nodes_ds, nodes):
         if node_ids.size > 0:
             missing_ids = sorted(set(nodes_ds[:]) - set(node_ids))
             if missing_ids:
-                errors.append(
-                    BluepySnapValidationError.fatal(
-                        f"{nodes_ds.name} misses node ids in its node population: {missing_ids}"
+                with _np_printoptions():
+                    errors.append(
+                        BluepySnapValidationError.fatal(
+                            f"{nodes_ds.name} misses node ids in its node population: {missing_ids}"
+                        )
                     )
-                )
         elif f"nodes/{node_population_name}" in h5f:
             errors.append(
                 BluepySnapValidationError.fatal(
@@ -386,12 +395,13 @@ def _check(indices, nodes_ds):
                 edges_range = node_to_edges_ranges[nodes_range[0] : nodes_range[1]][0]
                 edge_node_ids = list(set(nodes_ds[edges_range[0] : edges_range[1]]))
                 if len(edge_node_ids) > 1 or edge_node_ids[0] != node_id:
-                    errors.append(
-                        BluepySnapValidationError.fatal(
-                            f"Population {population.file.filename} edges {edge_node_ids} have "
-                            f"node ids {edges_range} instead of single id {node_id}"
+                    with _np_printoptions():
+                        errors.append(
+                            BluepySnapValidationError.fatal(
+                                f"Population {population.file.filename} edges {edge_node_ids} have "
+                                f"node ids {edges_range} instead of single id {node_id}"
+                            )
                         )
-                    )
 
     errors = []
     indices = population["indices"]
diff --git a/bluepysnap/edges/edge_population.py b/bluepysnap/edges/edge_population.py
index bbdeb459..114b39de 100644
--- a/bluepysnap/edges/edge_population.py
+++ b/bluepysnap/edges/edge_population.py
@@ -506,7 +506,7 @@ def _optimal_direction():
             # np.stack(uint64, int64) -> float64
             connected_node_ids_with_count = connected_node_ids_with_count.astype(np.uint32)
             if secondary_node_ids is not None:
-                mask = np.in1d(
+                mask = np.isin(
                     connected_node_ids_with_count[:, 0], secondary_node_ids, assume_unique=True
                 )
                 connected_node_ids_with_count = connected_node_ids_with_count[mask]
diff --git a/tests/test_edges.py b/tests/test_edges.py
index 7cab4001..114e5899 100644
--- a/tests/test_edges.py
+++ b/tests/test_edges.py
@@ -233,9 +233,9 @@ def test_get(self):
         tested = self.test_obj.get(ids, properties=["other2", "other1", "@source_node"])
         expected = pd.DataFrame(
             {
-                "other2": np.array([np.NaN, np.NaN, np.NaN, np.NaN, 10, 11, 12, 13], dtype=float),
+                "other2": np.array([np.nan, np.nan, np.nan, np.nan, 10, 11, 12, 13], dtype=float),
                 "other1": np.array(
-                    [np.NaN, np.NaN, np.NaN, np.NaN, "A", "B", "C", "D"], dtype=object
+                    [np.nan, np.nan, np.nan, np.nan, "A", "B", "C", "D"], dtype=object
                 ),
                 "@source_node": np.array([2, 0, 0, 2, 2, 0, 0, 2], dtype=int),
             },
@@ -329,7 +329,7 @@ def test_get(self):
         tested = self.test_obj.get(ids, properties="other2")
         expected = pd.DataFrame(
             {
-                "other2": np.array([np.NaN, np.NaN, np.NaN, np.NaN, 10, 11, 12, 13], dtype=float),
+                "other2": np.array([np.nan, np.nan, np.nan, np.nan, 10, 11, 12, 13], dtype=float),
             },
             index=pd.MultiIndex.from_tuples(
                 [
diff --git a/tests/test_nodes.py b/tests/test_nodes.py
index c135bfe9..d3f4da7e 100644
--- a/tests/test_nodes.py
+++ b/tests/test_nodes.py
@@ -286,8 +286,8 @@ def test_get(self):
         tested = pd.concat(df for _, df in tested)
         expected = pd.DataFrame(
             {
-                "other2": np.array([np.NaN, np.NaN, np.NaN, 10, 11, 12, 13], dtype=float),
-                "other1": np.array([np.NaN, np.NaN, np.NaN, "A", "B", "C", "D"], dtype=object),
+                "other2": np.array([np.nan, np.nan, np.nan, 10, 11, 12, 13], dtype=float),
+                "other1": np.array([np.nan, np.nan, np.nan, "A", "B", "C", "D"], dtype=object),
                 "layer": np.array([2, 6, 6, 7, 8, 8, 2], dtype=int),
             },
             index=pd.MultiIndex.from_tuples(
@@ -373,7 +373,7 @@ def test_get(self):
         tested = self.test_obj.get(properties="other2")
         expected = pd.DataFrame(
             {
-                "other2": np.array([np.NaN, np.NaN, np.NaN, 10, 11, 12, 13], dtype=float),
+                "other2": np.array([np.nan, np.nan, np.nan, 10, 11, 12, 13], dtype=float),
             },
             index=pd.MultiIndex.from_tuples(
                 [