diff --git a/src/gui/game/GameView.cpp b/src/gui/game/GameView.cpp
index 48f2c91d..d3b7b62c 100644
--- a/src/gui/game/GameView.cpp
+++ b/src/gui/game/GameView.cpp
@@ -2544,6 +2544,10 @@ void GameView::OnDraw()
        //FPS and some version info
        StringBuilder fpsInfo;
        fpsInfo << Format::Precision(2) << "FPS: " << ui::Engine::Ref().GetFps();
+       if (sim)
+       {
+           fpsInfo << " Threads: " << sim->GetActiveSimulationThreads() << "/" << sim->GetThreadCount();
+       }

        if (showDebug)
        {
diff --git a/src/simulation/Simulation.cpp b/src/simulation/Simulation.cpp
index 89b169a3..d39efe32 100644
--- a/src/simulation/Simulation.cpp
+++ b/src/simulation/Simulation.cpp
@@ -21,9 +21,44 @@
 #include <numbers>
 #include <set>
 #include <stack>
+#include <mutex>
+#include <thread>
+#include <condition_variable>
+#include <atomic>

 namespace
 {
+   struct SimulationImpl;
+
+   struct SimulationThreadPool
+   {
+       struct Worker
+       {
+           std::thread thread;
+           std::mutex mutex;
+           std::condition_variable cv;
+           bool hasWork = false;
+           bool stop = false;
+
+           SimulationImpl *sim = nullptr;
+           int start = 0;
+           int end = 0;
+       };
+
+       std::vector<std::unique_ptr<Worker>> workers;
+       std::mutex joinMutex;
+       std::condition_variable joinCv;
+       std::atomic<int> activeWorkers{0};
+
+   public:
+       SimulationThreadPool();
+       ~SimulationThreadPool();
+       void RunParallel(SimulationImpl *sim, int start, int end, int active);
+       int GetThreadCount() const { return workers.size(); }
+   };
+
+   static SimulationThreadPool globalThreadPool;
+
    struct SimulationImpl : public Simulation
    {
        struct Neighbourhood
@@ -39,7 +74,117 @@ namespace
        bool TransitionPhase(int i, const Neighbourhood &neighbourhood);

        void UpdateParticles(int start, int end) final override;
+       void UpdateParticlesSerial(int start, int end);
+       int GetThreadCount() const final override;
+       int GetActiveSimulationThreads() const final override;
    };
+
+   SimulationThreadPool::SimulationThreadPool()
+   {
+       int numThreads = std::thread::hardware_concurrency();
+       if (numThreads < 1)
+           numThreads = 1;
+       for (int i = 0; i < numThreads; ++i)
+       {
+           auto w = std::make_unique<Worker>();
+           w->thread = std::thread([this, workerPtr = w.get()]() {
+               while (true)
+               {
+                   std::unique_lock<std::mutex> lock(workerPtr->mutex);
+                   workerPtr->cv.wait(lock, [workerPtr]() { return workerPtr->hasWork || workerPtr->stop; });
+                   if (workerPtr->stop)
+                   {
+                       break;
+                   }
+
+                   // Do work
+                   workerPtr->sim->UpdateParticlesSerial(workerPtr->start, workerPtr->end);
+
+                   workerPtr->hasWork = false;
+                   lock.unlock();
+
+                   if (--activeWorkers == 0)
+                   {
+                       joinCv.notify_one();
+                   }
+               }
+           });
+           workers.push_back(std::move(w));
+       }
+   }
+
+   SimulationThreadPool::~SimulationThreadPool()
+   {
+       for (auto &w : workers)
+       {
+           {
+               std::lock_guard<std::mutex> lock(w->mutex);
+               w->stop = true;
+           }
+           w->cv.notify_one();
+           if (w->thread.joinable())
+           {
+               w->thread.join();
+           }
+       }
+   }
+
+   void SimulationThreadPool::RunParallel(SimulationImpl *sim, int start, int end, int active)
+   {
+       int totalRange = std::min(end, active) - start;
+       if (totalRange <= 0)
+       {
+           return;
+       }
+
+       int numWorkers = workers.size();
+       int chunkSize = (totalRange + numWorkers - 1) / numWorkers;
+       if (chunkSize < 1000)
+       {
+           // Too small, run sequentially on main thread
+           sim->UpdateParticlesSerial(start, end);
+           return;
+       }
+
+       activeWorkers = numWorkers;
+
+       for (int i = 0; i < numWorkers; ++i)
+       {
+           auto &w = workers[i];
+           int wStart = start + i * chunkSize;
+           int wEnd = std::min(wStart + chunkSize, end);
+
+           {
+               std::lock_guard<std::mutex> lock(w->mutex);
+               w->sim = sim;
+               w->start = wStart;
+               w->end = wEnd;
+               w->hasWork = true;
+           }
+           w->cv.notify_one();
+       }
+
+       // Wait for all workers to finish
+       std::unique_lock<std::mutex> lock(joinMutex);
+       joinCv.wait(lock, [this]() { return activeWorkers == 0; });
+   }
+
+   int SimulationImpl::GetThreadCount() const
+   {
+       return globalThreadPool.GetThreadCount();
+   }
+
+   int SimulationImpl::GetActiveSimulationThreads() const
+   {
+       if (ensureDeterminism)
+           return 1;
+       int totalRange = parts.active;
+       int numWorkers = globalThreadPool.GetThreadCount();
+       int chunkSize = (totalRange + numWorkers - 1) / numWorkers;
+       if (chunkSize < 1000)
+           return 1;
+       return numWorkers;
+   }
 }

 static float remainder_p(float x, float y)
@@ -1245,6 +1390,7 @@ int Simulation::try_move(int i, int x, int y, int nx, int ny)
            parts[ID(r)].tmp = (int)((parts[ID(r)].temp-73.15f)/100+1);
            if (parts[ID(r)].tmp>=CHANNELS) parts[ID(r)].tmp = CHANNELS-1;
            else if (parts[ID(r)].tmp<0) parts[ID(r)].tmp = 0;
+           std::lock_guard<std::recursive_mutex> lock(simulationMutex);
            for ( nnx=0; nnx<80; nnx++)
                if (!portalp[parts[ID(r)].tmp][count][nnx].type)
                {
@@ -1754,6 +1900,7 @@ Simulation::GetNormalResult Simulation::get_normal_interp<false, const Simulatio

 void Simulation::kill_part(int i)//kills particle number i
 {
+   std::lock_guard<std::recursive_mutex> lock(simulationMutex);
    if (i < 0 || i >= NPART)
        return;

@@ -1797,6 +1944,7 @@ void Parts::Free(int i)
 // Returns true if the particle was killed
 bool Simulation::part_change_type(int i, int x, int y, int t)
 {
+   std::lock_guard<std::recursive_mutex> lock(simulationMutex);
    if (x<0 || y<0 || x>=XRES || y>=YRES || i>=NPART || t<0 || t>=PT_NUM || !parts[i].type)
        return false;

@@ -1842,6 +1990,7 @@ bool Simulation::part_change_type(int i, int x, int y, int t)
 //tv = Type (PMAPBITS bits) + Var (32-PMAPBITS bits), var is usually 0
 int Simulation::create_part(int p, int x, int y, int t, int v)
 {
+   std::lock_guard<std::recursive_mutex> lock(simulationMutex);
    int i, oldType = PT_NONE;

    auto &sd = SimulationData::CRef();
@@ -2290,6 +2439,18 @@ SimulationImpl::Neighbourhood SimulationImpl::GetNeighbourhood(int i) const
 }

 void SimulationImpl::UpdateParticles(int start, int end)
+{
+   if (ensureDeterminism)
+   {
+       UpdateParticlesSerial(start, end);
+   }
+   else
+   {
+       globalThreadPool.RunParallel(this, start, end, parts.active);
+   }
+}
+
+void SimulationImpl::UpdateParticlesSerial(int start, int end)
 {
    //the main particle loop function, goes over all particles.
    auto &sd = SimulationData::CRef();
diff --git a/src/simulation/Simulation.h b/src/simulation/Simulation.h
index b401999c..fcb0e32a 100644
--- a/src/simulation/Simulation.h
+++ b/src/simulation/Simulation.h
@@ -20,6 +20,7 @@
 #include <array>
 #include <memory>
 #include <optional>
+#include <mutex>

 constexpr int CHANNELS = int(MAX_TEMP - 73) / 100 + 2;

@@ -269,6 +270,8 @@ public:
    virtual ~Simulation();

    void EnableNewtonianGravity(bool enable);
+   virtual int GetThreadCount() const { return 1; }
+   virtual int GetActiveSimulationThreads() const { return 1; }

    FrameTime *frameTime = nullptr;

@@ -280,4 +283,7 @@ private:
    void ResetNewtonianGravity(GravityInput newGravIn, GravityOutput newGravOut);
    void DispatchNewtonianGravity();
    void UpdateGravityMask();
+
+public:
+   std::recursive_mutex simulationMutex;
 };