diff --git a/avida-core/source/actions/PopulationActions.cc b/avida-core/source/actions/PopulationActions.cc
index 1b9ce83f0..947c0f259 100644
--- a/avida-core/source/actions/PopulationActions.cc
+++ b/avida-core/source/actions/PopulationActions.cc
@@ -57,6 +57,7 @@
 #include <limits>
 #include <numeric>
 #include <set>
+#include <vector>
 
 #include "stdlib.h"
 
@@ -5270,49 +5271,65 @@ class cActionKillDemesHighestParasiteLoad : public cAction
   {
     cPopulation& pop = m_world->GetPopulation();
     const int num_demes = pop.GetNumDemes();
+    std::vector<int> deme_indices(num_demes);
+    std::iota(std::begin(deme_indices), std::end(deme_indices), 0);
+
+    struct HasAny {
+      cPopulation& pop;
+      HasAny(cPopulation& pop) : pop(pop) {}
+      bool operator()(const int d){ return not pop.GetDeme(d).IsEmpty(); }
+    };
+    const int num_eligible = std::count_if(
+      std::begin(deme_indices), std::end(deme_indices),
+      HasAny(pop)
+    );
+    const int binomial_draw = ctx.GetRandom().GetRandBinomial(
+      num_eligible,
+      m_killprob
+    );
+    const int kill_quota = std::min(m_killmax, binomial_draw);
+    if (kill_quota != binomial_draw) {
+      std::cout << "warning: capped kill quota at " << kill_quota << " from " << binomial_draw << " binomial sample with " << num_eligible << " eligible and kill prob " << m_killprob << std::endl;
+    }
+
+    struct GetParasiteLoad {
+      cPopulation& pop;
+      GetParasiteLoad(cPopulation& pop) : pop(pop) {}
+      double operator()(const int d){ return pop.GetDeme(d).GetParasiteLoad(); }
+    };
     std::vector<double> parasite_loads(num_demes);
-    int num_eligible = 0;
+    std::transform(
+      std::begin(deme_indices), std::end(deme_indices),
+      std::begin(parasite_loads),
+      GetParasiteLoad(pop)
+    );
 
-    const int deme_size = m_world->GetConfig().WORLD_X.Get() * (m_world->GetConfig().WORLD_Y.Get() / num_demes);
-    const double smudge_delta = 0.09 / deme_size;
-    int smudge_index = ctx.GetRandom().GetInt(0, num_demes - 1);
-    for (int d = 0; d < num_demes; d++)
-    {
-      cDeme &deme = pop.GetDeme(d);
-      if (not deme.IsEmpty())
-      {
-        num_eligible++;
-        const auto parasite_load = deme.GetParasiteLoad();
-        if (parasite_load == 0.0) continue;
-        // need to guarantee that parasite_loads are distinct to set threshold
-        parasite_loads[d] = parasite_load + smudge_delta * smudge_index;
-        ++smudge_index;
-        if (smudge_index >= num_demes) smudge_index -= num_demes;
+    struct Comp {
+      std::vector<double>& loads;
+      Comp(std::vector<double> &loads) : loads(loads) {}
+      bool operator()(const int d1, const int d2) {
+        return loads[d1] > loads[d2];
       }
-    }
-
-    const int binomial_draw = ctx.GetRandom().GetRandBinomial(
-        num_eligible,
-        m_killprob
+    };
+    std::partial_sort(
+        std::begin(deme_indices),
+        std::next(std::begin(deme_indices), kill_quota),
+        std::end(deme_indices),
+        Comp(parasite_loads)
     );
-    const int kill_quota = std::min(binomial_draw, m_killmax);
-    if (kill_quota == 0) return;
-
-    std::vector<double> top_n(kill_quota);
-    const auto partial_sort_end = std::partial_sort_copy(
-      parasite_loads.begin(), parasite_loads.end(),
-      top_n.begin(), top_n.end(),
-      std::greater<int>()
+
+    struct DoKill {
+      cPopulation& pop;
+      cAvidaContext& ctx;
+      DoKill(cPopulation& pop, cAvidaContext& ctx) : pop(pop), ctx(ctx) {}
+      void operator()(const int d) { pop.GetDeme(d).KillAll(ctx); }
+    };
+    std::for_each(
+      std::begin(deme_indices),
+      std::next(std::begin(deme_indices), kill_quota),
+      DoKill(pop, ctx)
     );
-    const auto kill_thresh = *std::prev(partial_sort_end);
-    for (int d = 0; d < num_demes; d++)
-    {
-      if (parasite_loads[d] and parasite_loads[d] >= kill_thresh)
-      {
-        cDeme &deme = pop.GetDeme(d);
-        deme.KillAll(ctx);
-      }
-    }
+
 } // End Process()
 };
 
diff --git a/avida-core/source/main/cDeme.h b/avida-core/source/main/cDeme.h
index e557f638a..d7645f30d 100644
--- a/avida-core/source/main/cDeme.h
+++ b/avida-core/source/main/cDeme.h
@@ -187,7 +187,9 @@ class cDeme
 
   int GetNumParasites() const;
   double GetParasiteLoad() const {
-    return static_cast<double>(GetNumParasites()) / GetOrgCount();
+    const auto org_count = GetOrgCount();
+    if (org_count == 0) return 0.0;
+    else return static_cast<double>(GetNumParasites()) / GetOrgCount();
   }
 
   void UpdateParasiteMemoryScore(const double decay) {
diff --git a/avida-core/source/main/cPopulation.cc b/avida-core/source/main/cPopulation.cc
index 592b83377..2dfbdd54d 100644
--- a/avida-core/source/main/cPopulation.cc
+++ b/avida-core/source/main/cPopulation.cc
@@ -1211,6 +1211,12 @@ bool cPopulation::ActivateParasite(cOrganism* host, Systematics::UnitPtr parent,
   ){
     cDeme& deme = GetDeme(m_world->GetMigrationMatrix().GetProbabilisticDemeID(host_cell.GetDemeID(), m_world->GetRandom(),true));
     const int infection_mode = m_world->GetConfig().DEMES_PARASITE_MIGRATION_TARGET_SELECTION_METHOD.Get();
+    const auto immune_thresh = m_world->GetConfig().DEMES_PARASITE_MIGRATION_MEMORY_SCORE_PROTECTIVE_THRESHOLD.Get();
+    if (
+      immune_thresh != 0.0
+      && deme.GetParasiteMemoryScore() >= immune_thresh
+    ) return false;
+
     if (infection_mode == 0) {
     // Implementation #1 - Picks randomly of ALL cells in to-deme and then finds if the one it chose was occupied
     // -- Not ensured to infect an individual