diff --git a/tests/test_other_functions.py b/tests/test_other_functions.py
index 0978673..82a1e2d 100644
--- a/tests/test_other_functions.py
+++ b/tests/test_other_functions.py
@@ -161,7 +161,7 @@ def test_scenario_write_and_read():
     W.adddemand_point2point(0.5, 0.5, 2.5, 2.5, 0, 1800, volume=100)
     W.adddemand_area2area(0.5, 0.5, 2, 2.5, 2.5, 2, 0, 1800, volume=100)
 
-    W.save_scenario("out/test_grid.pkl")
+    W.save_scenario("out/test_grid.uxsim_scenario")
 
     W.exec_simulation()
     W.analyzer.print_simple_stats()
@@ -181,7 +181,7 @@ def test_scenario_write_and_read():
         random_seed=42
     )
 
-    W2.load_scenario("out/test_grid.pkl")
+    W2.load_scenario("out/test_grid.uxsim_scenario")
 
     W2.exec_simulation()
 
@@ -192,6 +192,79 @@ def test_scenario_write_and_read():
     
     assert df1["total_travel_time"][0] == df2["total_travel_time"][0]
 
+def test_scenario_write_and_read_areas():
+    ##################################
+    # Iter 1
+    ##################################
+
+    W = World(
+        name="",
+        deltan=10,
+        tmax=4000,
+        print_mode=1, save_mode=1, show_mode=0,
+        random_seed=42,
+    )
+
+    n_nodes = 5
+    imax = n_nodes
+    jmax = n_nodes
+    nodes = {}
+    for i in range(imax):
+        for j in range(jmax):
+            nodes[i,j] = W.addNode(f"n{(i,j)}", i, j, flow_capacity=1.6)
+
+    links = {}
+    for i in range(imax):
+        for j in range(jmax):
+            if i != imax-1:
+                links[i,j,i+1,j] = W.addLink(f"l{(i,j,i+1,j)}", nodes[i,j], nodes[i+1,j], length=1000)
+            if i != 0:
+                links[i,j,i-1,j] = W.addLink(f"l{(i,j,i-1,j)}", nodes[i,j], nodes[i-1,j], length=1000)
+            if j != jmax-1:
+                links[i,j,i,j+1] = W.addLink(f"l{(i,j,i,j+1)}", nodes[i,j], nodes[i,j+1], length=1000)
+            if j != 0:
+                links[i,j,i,j-1] = W.addLink(f"l{(i,j,i,j-1)}", nodes[i,j], nodes[i,j-1], length=1000)
+
+
+    areas = {
+        "areaN": [nodes[0,0], nodes[0, n_nodes-1]],
+        "areaS": [nodes[n_nodes-1,0], nodes[n_nodes-1, n_nodes-1]],
+        "areaNW": [nodes[0,0]],
+        "areaSE": [nodes[n_nodes-1, n_nodes-1]]
+    }
+
+    W.adddemand_areas2areas(areas["areaN"], areas["areaS"], 0, 3000, volume=7000)
+
+    W.save_scenario("out/test_area.uxsim_scenario")
+
+    W.exec_simulation()
+    W.analyzer.print_simple_stats()
+
+    df1 = W.analyzer.basic_to_pandas()
+    print(df1)
+
+    ##################################
+    # Iter 2
+    ##################################
+
+    W2 = World(
+        name="",
+        deltan=10,
+        tmax=4000,
+        print_mode=1, save_mode=1, show_mode=0,
+        random_seed=42,
+    )
+
+    W2.load_scenario("out/test_area.uxsim_scenario")
+
+    W2.exec_simulation()
+    W2.analyzer.print_simple_stats()
+
+    df2 = W2.analyzer.basic_to_pandas()
+    print(df2)
+
+    assert df1["total_travel_time"][0] == df2["total_travel_time"][0]
+
 def test_k_shortest_path():
     W = World(
         name="",    # Scenario name
@@ -420,7 +493,7 @@ def test_area():
         "areaSE": [nodes[n_nodes-1, n_nodes-1]]
     }
 
-    W.adddemand_nodes2nodes(areas["areaN"], areas["areaS"], 0, 3000, volume=7000)
+    W.adddemand_areas2areas(areas["areaN"], areas["areaS"], 0, 3000, volume=7000)
 
     W.exec_simulation()
     W.analyzer.print_simple_stats()
@@ -436,4 +509,4 @@ def test_area():
     assert df["average_shortest_distance"][(df["origin_area"] == "areaN") & (df["destination_area"] == "areaS")].values[0] == 6000.0
     assert df["total_trips"][(df["origin_area"] == "areaNW") & (df["destination_area"] == "areaSE")].values[0] == 1750
     assert df["average_free_travel_time"][(df["origin_area"] == "areaNW") & (df["destination_area"] == "areaSE")].values[0] == 400.0
-    assert df["average_shortest_distance"][(df["origin_area"] == "areaNW") & (df["destination_area"] == "areaSE")].values[0] == 8000.0
\ No newline at end of file
+    assert df["average_shortest_distance"][(df["origin_area"] == "areaNW") & (df["destination_area"] == "areaSE")].values[0] == 8000.0
diff --git a/uxsim/uxsim.py b/uxsim/uxsim.py
index 7fc8a24..5541f70 100644
--- a/uxsim/uxsim.py
+++ b/uxsim/uxsim.py
@@ -1802,10 +1802,9 @@ def adddemand_area2area(W, x_orig, y_orig,  radious_orig, x_dest, y_dest, radiou
             for d in dests:
                 W.adddemand(o, d, t_start, t_end, flow, volume, attribute, direct_call=False)
     
-    #@demand_info_record: TODO: activate this
-    def adddemand_nodes2nodes(W, origs, dests, t_start, t_end, flow=-1, volume=-1, attribute=None, direct_call=True):
+    def adddemand_areas2areas(W, origs, dests, t_start, t_end, flow=-1, volume=-1, attribute=None):
         """
-        Generate vehicles by specifying time-dependent origin-destination demand by specifying sets of nodes.
+        Generate vehicles by specifying time-dependent origin-destination demand by specifying areas (i.e., sets of nodes).
 
         Parameters
         ----------
@@ -1844,7 +1843,7 @@ def adddemand_nodes2nodes(W, origs, dests, t_start, t_end, flow=-1, volume=-1, a
             volume = volume/(len(origs)*len(dests))
         for o in origs:
             for d in dests:
-                W.adddemand(o, d, t_start, t_end, flow, volume, attribute, direct_call=False)
+                W.adddemand(o, d, t_start, t_end, flow, volume, attribute, direct_call=True)
 
     def finalize_scenario(W, tmax=None):
         """