diff --git a/src/dsm/headers/Dynamics.hpp b/src/dsm/headers/Dynamics.hpp
index dc984616e..50c803b08 100644
--- a/src/dsm/headers/Dynamics.hpp
+++ b/src/dsm/headers/Dynamics.hpp
@@ -25,6 +25,8 @@
 #include <filesystem>
 #include <functional>
 
+#include <tbb/tbb.h>
+
 #include "Network.hpp"
 #include "../utility/TypeTraits/is_agent.hpp"
 #include "../utility/TypeTraits/is_itinerary.hpp"
@@ -67,6 +69,7 @@ namespace dsm {
     network_t m_graph;
 
   protected:
+    tbb::task_arena m_taskArena;
     Time m_time;
     std::mt19937_64 m_generator;
 
@@ -96,5 +99,6 @@ namespace dsm {
     if (seed.has_value()) {
       m_generator.seed(*seed);
     }
+    m_taskArena.initialize();
   }
 };  // namespace dsm
diff --git a/src/dsm/headers/RoadDynamics.hpp b/src/dsm/headers/RoadDynamics.hpp
index ca2efb38f..6be4e699c 100644
--- a/src/dsm/headers/RoadDynamics.hpp
+++ b/src/dsm/headers/RoadDynamics.hpp
@@ -52,8 +52,8 @@ namespace dsm {
     std::unordered_map<Id, std::array<unsigned long long, 4>> m_turnCounts;
     std::unordered_map<Id, std::array<long, 4>> m_turnMapping;
     std::unordered_map<Id, double> m_streetTails;
-    std::vector<std::pair<double, double>> m_travelDTs;
-    std::vector<Id> m_agentsToRemove;
+    tbb::concurrent_vector<std::pair<double, double>> m_travelDTs;
+    tbb::concurrent_vector<Id> m_agentsToRemove;
     Time m_previousOptimizationTime, m_previousSpireTime;
 
   private:
@@ -645,16 +645,16 @@ namespace dsm {
       auto const& nextStreet{
           this->graph().edge(this->agents().at(agentId)->nextStreetId().value())};
       if (!(nextStreet->isFull())) {
-        if (this->agents().at(agentId)->streetId().has_value()) {
-          const auto streetId = this->agents().at(agentId)->streetId().value();
-          auto delta = nextStreet->angle() - this->graph().edge(streetId)->angle();
-          if (delta > std::numbers::pi) {
-            delta -= 2 * std::numbers::pi;
-          } else if (delta < -std::numbers::pi) {
-            delta += 2 * std::numbers::pi;
-          }
-          m_increaseTurnCounts(streetId, delta);
-        }
+        // if (this->agents().at(agentId)->streetId().has_value()) {
+        //   const auto streetId = this->agents().at(agentId)->streetId().value();
+        //   auto delta = nextStreet->angle() - this->graph().edge(streetId)->angle();
+        //   if (delta > std::numbers::pi) {
+        //     delta -= 2 * std::numbers::pi;
+        //   } else if (delta < -std::numbers::pi) {
+        //     delta += 2 * std::numbers::pi;
+        //   }
+        //   m_increaseTurnCounts(streetId, delta);
+        // }
         roundabout.dequeue();
         this->agents().at(agentId)->setStreetId(nextStreet->id());
         this->setAgentSpeed(agentId);
@@ -1078,45 +1078,64 @@ namespace dsm {
   void RoadDynamics<delay_t>::evolve(bool reinsert_agents) {
     // move the first agent of each street queue, if possible, putting it in the next node
     bool const bUpdateData =
-        m_dataUpdatePeriod.has_value() && this->time() % m_dataUpdatePeriod.value() == 0;
-    auto const N{this->graph().nNodes()};
-    std::for_each(this->graph().nodes().cbegin(),
-                  this->graph().nodes().cend(),
-                  [&](const auto& pair) {
-                    for (auto const& sourceId :
-                         this->graph().adjacencyMatrix().getCol(pair.first)) {
-                      auto const streetId = sourceId * N + pair.first;
-                      auto const& pStreet{this->graph().edge(streetId)};
-                      if (bUpdateData) {
-                        m_streetTails[streetId] += pStreet->nExitingAgents();
-                      }
-                      for (auto i = 0; i < pStreet->transportCapacity(); ++i) {
-                        m_evolveStreet(pStreet, reinsert_agents);
-                      }
-                    }
-                  });
+        m_dataUpdatePeriod.has_value() && this->m_time % m_dataUpdatePeriod.value() == 0;
+    auto const numNodes{this->graph().nNodes()};
+    const auto& nodes =
+        this->graph().nodes();  // assuming a container with contiguous indices
+    const unsigned int concurrency = std::thread::hardware_concurrency();
+    // Calculate a grain size to partition the nodes into roughly "concurrency" blocks
+    const size_t grainSize = std::max(size_t(1), numNodes / concurrency);
+    this->m_taskArena.execute([&] {
+      tbb::parallel_for(tbb::blocked_range<size_t>(0, numNodes, grainSize),
+                        [&](const tbb::blocked_range<size_t>& range) {
+                          for (size_t i = range.begin(); i != range.end(); ++i) {
+                            const auto& pNode = nodes.at(i);
+                            for (auto const& sourceId :
+                                 this->graph().adjacencyMatrix().getCol(pNode->id())) {
+                              auto const streetId = sourceId * numNodes + pNode->id();
+                              auto const& pStreet = this->graph().edge(streetId);
+
+                              if (bUpdateData) {
+                                m_streetTails[streetId] += pStreet->nExitingAgents();
+                              }
+
+                              for (auto j = 0; j < pStreet->transportCapacity(); ++j) {
+                                this->m_evolveStreet(pStreet, reinsert_agents);
+                              }
+                            }
+                          }
+                        });
+    });
     std::for_each(this->m_agentsToRemove.cbegin(),
                   this->m_agentsToRemove.cend(),
                   [this](const auto& agentId) { this->removeAgent(agentId); });
     m_agentsToRemove.clear();
     // Move transport capacity agents from each node
-    std::for_each(this->graph().nodes().cbegin(),
-                  this->graph().nodes().cend(),
-                  [&](const auto& pair) {
-                    for (auto i = 0; i < pair.second->transportCapacity(); ++i) {
-                      m_evolveNode(pair.second);
-                    }
-                    if (pair.second->isTrafficLight()) {
-                      auto& tl = dynamic_cast<TrafficLight&>(*pair.second);
-                      ++tl;
-                    }
-                  });
+    this->m_taskArena.execute([&] {
+      tbb::parallel_for(tbb::blocked_range<size_t>(0, numNodes, grainSize),
+                        [&](const tbb::blocked_range<size_t>& range) {
+                          for (size_t i = range.begin(); i != range.end(); ++i) {
+                            const auto& pNode = nodes.at(i);
+
+                            // Evolve the node based on its transport capacity.
+                            for (auto j = 0; j < pNode->transportCapacity(); ++j) {
+                              this->m_evolveNode(pNode);
+                            }
+
+                            // If the node is a traffic light, increment it.
+                            if (pNode->isTrafficLight()) {
+                              auto& tl = dynamic_cast<TrafficLight&>(*pNode);
+                              ++tl;
+                            }
+                          }
+                        });
+    });
     // cycle over agents and update their times
-    std::for_each(this->agents().cbegin(),
-                  this->agents().cend(),
-                  [this](auto const& pair) { m_evolveAgent(pair.second); });
-
-    Dynamics<RoadNetwork>::m_evolve();
+    std::for_each(this->agents().begin(), this->agents().end(), [&](auto const& pair) {
+      this->m_evolveAgent(pair.second);
+    });
+    // increment time simulation
+    ++this->m_time;
   }
 
   template <typename delay_t>