diff --git a/src/ivoc/ivocrand.cpp b/src/ivoc/ivocrand.cpp
index f06e899278..d3501913db 100644
--- a/src/ivoc/ivocrand.cpp
+++ b/src/ivoc/ivocrand.cpp
@@ -10,6 +10,7 @@
 #include "classreg.h"
 #include "oc2iv.h"
 #include "nrnisaac.h"
+#include "utils/enumerate.h"
 
 #include <vector>
 #include <ocnotify.h>
@@ -183,13 +184,11 @@ void RandomPlay::play() {
     *px_ = (*(r_->rand))();
 }
 void RandomPlay::list_remove() {
-    for (auto it = random_play_list_->begin(); it != random_play_list_->end(); ++it) {
-        if (*it == (RandomPlay*) this) {
-            // printf("RandomPlay %p removed from list cnt=%d i=%d %p\n", this, cnt, i);
-            random_play_list_->erase(it);
-            unref_deferred();
-            break;
-        }
+    if (auto it = std::find(random_play_list_->begin(), random_play_list_->end(), this);
+        it != random_play_list_->end()) {
+        // printf("RandomPlay %p removed from list cnt=%d i=%d %p\n", this, cnt, i);
+        random_play_list_->erase(it);
+        unref_deferred();
     }
 }
 void RandomPlay::update(Observable*) {
diff --git a/src/ivoc/ocpicker.cpp b/src/ivoc/ocpicker.cpp
index 76da6bee49..1271335d3d 100644
--- a/src/ivoc/ocpicker.cpp
+++ b/src/ivoc/ocpicker.cpp
@@ -1,6 +1,7 @@
 #include <../../nrnconf.h>
 #if HAVE_IV  // to end of file
 
+#include <algorithm>
 #include <OS/list.h>
 #include <InterViews/hit.h>
 #include <stdio.h>
@@ -103,16 +104,14 @@ void StandardPicker::event(const Event& e) {
 }
 
 void StandardPicker::unbind(int m, EventButton eb) {
-    long cnt = handlers_[m]->size();
-    for (long i = 0, j = 0; i < cnt; ++i) {
-        ButtonHandler* b = handlers_[m]->at(j);
+    for (auto& b: *handlers_[m]) {
         if (b->eb_ == Event::any || b->eb_ == eb) {
             delete b;
-            handlers_[m]->erase(handlers_[m]->begin() + j);
-        } else {
-            ++j;
+            b = nullptr;
         }
     }
+    handlers_[m]->erase(std::remove(handlers_[m]->begin(), handlers_[m]->end(), nullptr),
+                        handlers_[m]->end());
 }
 
 void StandardPicker::bind(int m, EventButton eb, OcHandler* h) {
diff --git a/src/ivoc/scene.cpp b/src/ivoc/scene.cpp
index 406846ad73..8011ae81eb 100644
--- a/src/ivoc/scene.cpp
+++ b/src/ivoc/scene.cpp
@@ -51,6 +51,7 @@
 #include "scenepic.h"
 #include "idraw.h"
 #include "ivoc.h"
+#include "utils/enumerate.h"
 
 #define Scene_Move_Text_   "MoveText Graph"
 #define Scene_ChangeColor_ "ChangeColor Graph"
@@ -311,10 +312,7 @@ Scene::~Scene() {
     // without first deleteing all the views.
     assert(views_->empty());
 
-    if (auto it = std::find(scene_list->begin(), scene_list->end(), this);
-        it != scene_list->end()) {
-        scene_list->erase(it);
-    }
+    erase_first(*scene_list, this);
     delete views_;
 }
 
@@ -335,14 +333,12 @@ void Scene::append_view(XYView* v) {
 }
 
 void Scene::remove_view(XYView* v) {
-    if (auto it = std::find(views_->begin(), views_->end(), v); it != views_->end()) {
-        views_->erase(it);
-    }
+    erase_first(*views_, v);
 }
 
 void Scene::dismiss() {
-    for (auto it = views_->rbegin(); it != views_->rend(); ++it) {
-        OcViewGlyph* g = (*it)->parent();
+    for (auto&& item: reverse(*views_)) {
+        OcViewGlyph* g = item->parent();
         if (g && g->has_window()) {
             g->window()->dismiss();
             g->window(NULL);
@@ -513,8 +509,8 @@ void Scene::replace(GlyphIndex index, Glyph* glyph) {
 }
 
 GlyphIndex Scene::glyph_index(const Glyph* g) {
-    for (std::size_t i = 0; i < info_->size(); ++i) {
-        if ((*info_)[i].glyph_ == g) {
+    for (const auto&& [i, info]: enumerate(*info_)) {
+        if (info.glyph_ == g) {
             return i;
         }
     }
@@ -536,7 +532,7 @@ void Scene::request(Requisition& req) const {
 
 void Scene::allocate(Canvas* c, const Allocation& a, Extension& ext) {
     // printf("Scene::allocate\n");
-    for (GlyphIndex index = 0; index < info_->size(); ++index) {
+    for (std::size_t index: range(*info_)) {
         check_allocation(index);
     }
     ext.set(c, a);
@@ -693,14 +689,12 @@ void Scene::pick(Canvas* c, const Allocation& a, int depth, Hit& h) {
     if (background_ != NULL) {
         background_->pick(c, a, depth, h);
     }
-    GlyphIndex count = info_->size();
     // pick with some extra epsilon in canvas coords
     Coord epsx = XYView::current_pick_view()->x_pick_epsilon();
     Coord epsy = XYView::current_pick_view()->y_pick_epsilon();
 
     bool are_fixed = false;
-    for (GlyphIndex index = 0; index < count; ++index) {
-        SceneInfo& info = (*info_)[index];
+    for (auto&& [index, info]: enumerate(*info_)) {
         if (info.status_ & SceneInfoFixed) {
             are_fixed = true;
         } else if (info.glyph_ != NULL && (info.status_ & SceneInfoShowing)) {
@@ -720,8 +714,7 @@ void Scene::pick(Canvas* c, const Allocation& a, int depth, Hit& h) {
         const Transformer& tv = XYView::current_pick_view()->s2o();
         float scx, scy, tmp;
         tv.matrix(scx, tmp, tmp, scy, tmp, tmp);
-        for (GlyphIndex index = 0; index < count; ++index) {
-            SceneInfo& info = (*info_)[index];
+        for (auto&& [index, info]: enumerate(*info_)) {
             if ((info.status_ & SceneInfoFixed) && info.glyph_ != NULL &&
                 (info.status_ & SceneInfoShowing)) {
                 Allocation a = info.allocation_;
@@ -752,13 +745,11 @@ void Scene::pick(Canvas* c, const Allocation& a, int depth, Hit& h) {
 }
 
 long Scene::scene_list_index(Scene* s) {
-    std::size_t cnt = scene_list->size();
-    for (std::size_t i = 0; i < cnt; ++i) {
-        if (s == scene_list->at(i)) {
-            return i;
-        }
+    auto it = std::find(scene_list->begin(), scene_list->end(), s);
+    if (it == scene_list->end()) {
+        return -1;
     }
-    return -1;
+    return std::distance(scene_list->begin(), it);
 }
 
 void Scene::save_all(std::ostream& o) {
diff --git a/src/ivoc/scenepic.cpp b/src/ivoc/scenepic.cpp
index 9b4b29d172..bd120ca726 100644
--- a/src/ivoc/scenepic.cpp
+++ b/src/ivoc/scenepic.cpp
@@ -16,6 +16,7 @@
 #include "apwindow.h"
 #include "utility.h"
 #include "oc2iv.h"
+#include "utils/enumerate.h"
 
 #define Scene_SceneMover_     "Translate Scene"
 #define Scene_SceneZoom_      "ZoomInOut Scene"
@@ -309,9 +310,7 @@ MenuItem* ScenePicker::add_radio_menu(const char* name, OcHandler* h, int tool,
 }
 
 long ScenePickerImpl::info_index(const char* name) {
-    std::size_t cnt = bil_->size();
-    for (std::size_t i = 0; i < cnt; ++i) {
-        ButtonItemInfo* b = bil_->at(i);
+    for (const auto&& [i, b]: enumerate(*bil_)) {
         if (strcmp(b->name_.string(), name) == 0) {
             return i;
         }
@@ -410,8 +409,8 @@ ScenePickerImpl::ScenePickerImpl(Scene* scene)
 ScenePickerImpl::~ScenePickerImpl() {
     Resource::unref(menu_);
     Resource::unref(tg_);
-    for (auto it = bil_->rbegin(); it != bil_->rend(); ++it) {
-        delete *it;
+    for (ButtonItemInfo* bii: *bil_) {
+        delete bii;
     }
     delete bil_;
 }
diff --git a/src/ivoc/symdir.cpp b/src/ivoc/symdir.cpp
index 46c55efbd9..2a06893c21 100644
--- a/src/ivoc/symdir.cpp
+++ b/src/ivoc/symdir.cpp
@@ -5,6 +5,7 @@
 #include <OS/string.h>
 #include <stdio.h>
 #include "ocobserv.h"
+#include "utils/enumerate.h"
 
 #include "nrniv_mf.h"
 #include "nrnoc2iv.h"
@@ -60,32 +61,16 @@ class SymDirectoryImpl: public Observer {
     void sort();
 };
 
-static int compare_entries(const void* k1, const void* k2) {
-    SymbolItem* e1 = *((SymbolItem**) k1);
-    SymbolItem* e2 = *((SymbolItem**) k2);
+static int compare_entries(const SymbolItem* e1, const SymbolItem* e2) {
     int i = strcmp(e1->name().string(), e2->name().string());
     if (i == 0) {
-        if (e1->array_index() > e2->array_index()) {
-            i = 1;
-        } else {
-            i = -1;
-        }
+        return e1->array_index() > e2->array_index();
     }
-    return i;
+    return i > 0;
 };
 
 void SymDirectoryImpl::sort() {
-    std::size_t cnt = symbol_lists_.size();
-    SymbolItem** slist = new SymbolItem*[cnt];
-    for (std::size_t i = 0; i < cnt; ++i) {
-        slist[i] = symbol_lists_[i];
-    }
-    qsort(slist, cnt, sizeof(SymbolItem*), compare_entries);
-    symbol_lists_.clear();
-    for (std::size_t i = 0; i < cnt; ++i) {
-        symbol_lists_.push_back(slist[i]);
-    }
-    delete[] slist;
+    std::sort(symbol_lists_.begin(), symbol_lists_.end(), compare_entries);
 }
 
 // SymDirectory
@@ -307,9 +292,8 @@ int SymDirectory::array_index(int i) const {
 }
 
 int SymDirectory::index(const String& name) const {
-    long cnt = count();
-    for (long i = 0; i < cnt; ++i) {
-        if (name == impl_->symbol_lists_[i]->name()) {
+    for (const auto&& [i, symbol]: enumerate(impl_->symbol_lists_)) {
+        if (name == symbol->name()) {
             return i;
         }
     }
@@ -451,8 +435,7 @@ void SymDirectoryImpl::load(int type) {
 }
 
 void SymDirectoryImpl::load(int type, Symlist* sl) {
-    Symbol* sym;
-    for (sym = sl->first; sym; sym = sym->next) {
+    for (Symbol* sym = sl->first; sym; sym = sym->next) {
         if (type == -1) {
             switch (sym->type) {
             case SECTION:
@@ -490,9 +473,8 @@ void SymDirectoryImpl::load_aliases() {
     IvocAliases* a = (IvocAliases*) obj_->aliases;
     if (!a)
         return;
-    for (const auto& kv: a->symtab_) {
-        Symbol* s = kv.second;
-        append(s, NULL, obj_);
+    for (const auto& [_, s]: a->symtab_) {
+        append(s, nullptr, obj_);
     }
 }
 
diff --git a/src/ivoc/xmenu.cpp b/src/ivoc/xmenu.cpp
index b8011df544..082db077ec 100644
--- a/src/ivoc/xmenu.cpp
+++ b/src/ivoc/xmenu.cpp
@@ -1,9 +1,9 @@
-
 #include <../../nrnconf.h>
 
 #include "oc2iv.h"
 #include "classreg.h"
 #include "gui-redirect.h"
+#include "utils/enumerate.h"
 
 #if HAVE_IV  // to end of file except for a few small fragments.
 
@@ -992,8 +992,8 @@ void HocPanel::map_window(int scroll) {
 // HocPanel
 static void var_freed(void* pd, int size) {
     if (hoc_panel_list) {
-        for (auto it = hoc_panel_list->rbegin(); it != hoc_panel_list->rend(); ++it) {
-            (*it)->check_valid_pointers(pd, size);
+        for (auto&& elem: reverse(*hoc_panel_list)) {
+            elem->check_valid_pointers(pd, size);
         }
     }
 }
@@ -1034,10 +1034,7 @@ HocPanel::~HocPanel() {
     for (auto& item: elist_) {
         item->HocItem::unref();
     }
-    if (auto it = std::find(hoc_panel_list->begin(), hoc_panel_list->end(), this);
-        it != hoc_panel_list->end()) {
-        hoc_panel_list->erase(it);
-    }
+    erase_first(*hoc_panel_list, this);
     ilist_.clear();
     ilist_.shrink_to_fit();
     elist_.clear();
@@ -1073,10 +1070,7 @@ void HocPanel::keep_updated(HocUpdateItem* hui, bool add) {
     if (add) {
         update_list_->push_back(hui);
     } else {
-        if (auto it = std::find(update_list_->begin(), update_list_->end(), hui);
-            it != update_list_->end()) {
-            update_list_->erase(it);
-        }
+        erase_first(*update_list_, hui);
     }
 }
 
@@ -2279,22 +2273,22 @@ void Oc::notifyHocValue() {
     ParseTopLevel ptl;
     ptl.save();
     if (hoc_panel_list) {
-        for (auto it = hoc_panel_list->rbegin(); it != hoc_panel_list->rend(); ++it) {
-            (*it)->notifyHocValue();
+        for (auto&& e: reverse(*hoc_panel_list)) {
+            e->notifyHocValue();
         }
     }
     ptl.restore();
 }
 
 void HocPanel::notifyHocValue() {
-    for (auto it = elist_.rbegin(); it != elist_.rend(); ++it) {
-        (*it)->update_hoc_item();
+    for (auto&& e: reverse(elist_)) {
+        e->update_hoc_item();
     }
 }
 
 void HocPanel::check_valid_pointers(void* v, int size) {
-    for (auto it = elist_.rbegin(); it != elist_.rend(); ++it) {
-        (*it)->check_pointer(v, size);
+    for (auto&& e: reverse(elist_)) {
+        e->check_pointer(v, size);
     }
 }
 
@@ -2321,8 +2315,8 @@ void HocVarLabel::check_pointer(void* v, int) {
 }
 
 void HocPanel::data_path(HocDataPaths* hdp, bool append) {
-    for (auto it = elist_.rbegin(); it != elist_.rend(); ++it) {
-        (*it)->data_path(hdp, append);
+    for (auto&& e: reverse(elist_)) {
+        e->data_path(hdp, append);
     }
 }
 
diff --git a/src/nrncvode/netcvode.cpp b/src/nrncvode/netcvode.cpp
index 53da0e10b4..2f9acf4fe3 100644
--- a/src/nrncvode/netcvode.cpp
+++ b/src/nrncvode/netcvode.cpp
@@ -22,6 +22,7 @@
 #include "shared/sundialsmath.h"
 #include "kssingle.h"
 #include "ocnotify.h"
+#include "utils/enumerate.h"
 #if HAVE_IV
 #include "ivoc.h"
 #include "glinerec.h"
@@ -3050,21 +3051,17 @@ static OcList* event_info_list_;  // netcon or point_process
 static void event_info_callback(const TQItem*, int);
 static void event_info_callback(const TQItem* q, int) {
     DiscreteEvent* d = (DiscreteEvent*) q->data_;
-    NetCon* nc;
-    PreSyn* ps;
-    SelfEvent* se;
-    int n = event_info_tvec_->size();
     switch (d->type()) {
     case NetConType:
         if (event_info_type_ == NetConType) {
-            nc = (NetCon*) d;
+            auto* nc = static_cast<NetCon*>(d);
             event_info_tvec_->push_back(q->t_);
             event_info_list_->append(nc->obj_);
         }
         break;
     case SelfEventType:
         if (event_info_type_ == SelfEventType) {
-            se = (SelfEvent*) d;
+            auto* se = static_cast<SelfEvent*>(d);
             event_info_tvec_->push_back(q->t_);
             event_info_flagvec_->push_back(se->flag_);
             event_info_list_->append(se->target_->ob);
@@ -3072,13 +3069,11 @@ static void event_info_callback(const TQItem* q, int) {
         break;
     case PreSynType:
         if (event_info_type_ == NetConType) {
-            ps = (PreSyn*) d;
-            for (auto it = ps->dil_.rbegin(); it != ps->dil_.rend(); ++it) {
-                nc = *it;
+            auto* ps = static_cast<PreSyn*>(d);
+            for (const auto& nc: reverse(ps->dil_)) {
                 double td = nc->delay_ - ps->delay_;
                 event_info_tvec_->push_back(q->t_ + td);
                 event_info_list_->append(nc->obj_);
-                ++n;
             }
         }
         break;
@@ -5774,9 +5769,8 @@ void nrnthread_get_trajectory_requests(int tid,
             int err = 0;
             if (pr->ith_ == tid) {
                 if (1) {  // buffered or per time step value return
-                    IvocVect* v = NULL;
                     if (pr->type() == TvecRecordType) {
-                        v = ((TvecRecord*) pr)->t_;
+                        IvocVect* v = ((TvecRecord*) pr)->t_;
                         err = trajec_buffered(nt,
                                               bsize,
                                               v,
@@ -5795,7 +5789,7 @@ void nrnthread_get_trajectory_requests(int tid,
                             n_trajec--;
                         }
                     } else if (pr->type() == YvecRecordType) {
-                        v = ((YvecRecord*) pr)->y_;
+                        IvocVect* v = ((YvecRecord*) pr)->y_;
                         err = trajec_buffered(nt,
                                               bsize,
                                               v,
@@ -5819,7 +5813,7 @@ void nrnthread_get_trajectory_requests(int tid,
                             if (bsize && !glr->v_) {
                                 glr->v_ = new IvocVect(bsize);
                             }
-                            v = glr->v_;
+                            IvocVect* v = glr->v_;
                             err = trajec_buffered(nt,
                                                   bsize,
                                                   v,
@@ -5837,14 +5831,11 @@ void nrnthread_get_trajectory_requests(int tid,
                                 n_trajec--;
                             }
                         } else {  // glr->gl_->name expression involves several range variables
-                            GLineRecordEData& ed = glr->pd_and_vec_;
                             int n = n_trajec;
-                            for (GLineRecordEData::iterator it = ed.begin(); it != ed.end(); ++it) {
-                                double* pd = (*it).first;
+                            for (auto&& [pd, v]: glr->pd_and_vec_) {
                                 assert(pd);
-                                v = (*it).second;
-                                if (bsize && v == NULL) {
-                                    v = (*it).second = new IvocVect(bsize);
+                                if (bsize && v == nullptr) {
+                                    v = new IvocVect(bsize);
                                 }
                                 // TODO avoid the conversion?
                                 err = trajec_buffered(nt,
@@ -6124,22 +6115,14 @@ void NetCvode::playrec_add(PlayRecord* pr) {  // called by PlayRecord constructo
 void NetCvode::playrec_remove(PlayRecord* pr) {  // called by PlayRecord destructor
                                                  // printf("NetCvode::playrec_remove %p\n", pr);
     playrec_change_cnt_ = 0;
-    if (auto it = std::find(prl_->begin(), prl_->end(), pr); it != prl_->end()) {
-        prl_->erase(it);
-    }
-    if (auto it = std::find(fixed_play_->begin(), fixed_play_->end(), pr);
-        it != fixed_play_->end()) {
-        fixed_play_->erase(it);
-    }
-    if (auto it = std::find(fixed_record_->begin(), fixed_record_->end(), pr);
-        it != fixed_record_->end()) {
-        fixed_record_->erase(it);
-    }
+    erase_first(*prl_, pr);
+    erase_first(*fixed_play_, pr);
+    erase_first(*fixed_record_, pr);
 }
 
 int NetCvode::playrec_item(PlayRecord* pr) {
-    for (std::size_t i = 0; i < prl_->size(); ++i) {
-        if ((*prl_)[i] == pr) {
+    for (const auto&& [i, e]: enumerate(*prl_)) {
+        if (e == pr) {
             return i;
         }
     }
diff --git a/src/nrniv/finithnd.cpp b/src/nrniv/finithnd.cpp
index 7ba8acd5c9..c9d237f935 100644
--- a/src/nrniv/finithnd.cpp
+++ b/src/nrniv/finithnd.cpp
@@ -19,6 +19,7 @@ Type 3 are at the very beginning of finitialize. ie structure changes
 #include <nrnoc2iv.h>
 #include <classreg.h>
 #include <objcmd.h>
+#include "utils/enumerate.h"
 
 class FInitialHandler {
   public:
@@ -106,10 +107,5 @@ FInitialHandler::FInitialHandler(int i, const char* s, Object* obj, Object* pyac
 
 FInitialHandler::~FInitialHandler() {
     delete stmt_;
-    for (auto it = fihlist_[type_].begin(); it != fihlist_[type_].end(); ++it) {
-        if ((*it) == this) {
-            fihlist_[type_].erase(it);
-            return;
-        }
-    }
+    erase_first(fihlist_[type_], this);
 }
diff --git a/src/nrniv/glinerec.cpp b/src/nrniv/glinerec.cpp
index 001550bd10..58cf63099e 100644
--- a/src/nrniv/glinerec.cpp
+++ b/src/nrniv/glinerec.cpp
@@ -18,6 +18,7 @@
 #include "vrecitem.h"
 #include "netcvode.h"
 #include "cvodeobj.h"
+#include "utils/enumerate.h"
 
 #if HAVE_IV  // to end of file
 #include "graph.h"
@@ -108,9 +109,9 @@ void GLineRecord::fill_pd1() {
 
 void GLineRecord::fill_pd() {
     // Get rid of old pd_and_vec_ info.
-    for (GLineRecordEData::iterator it = pd_and_vec_.begin(); it != pd_and_vec_.end(); ++it) {
-        if ((*it).second) {
-            delete (*it).second;
+    for (auto& [_, elem]: pd_and_vec_) {
+        if (elem) {
+            delete elem;
         }
     }
     pd_and_vec_.resize(0);
@@ -205,10 +206,8 @@ void GLineRecord::plot(int vecsz, double tstop) {
         ObjectContext obc(NULL);
         for (int i = 0; i < vecsz; ++i) {
             x->add(dt * i);
-            for (GLineRecordEData::iterator it = pd_and_vec_.begin(); it != pd_and_vec_.end();
-                 ++it) {
-                double* pd = (*it).first;
-                *pd = (*it).second->elem(i);
+            for (auto& [pd, elem]: pd_and_vec_) {
+                *pd = elem->elem(i);
             }
             gl_->plot();
         }
diff --git a/src/nrniv/nrndae.cpp b/src/nrniv/nrndae.cpp
index c383728011..4bbe44d8f6 100644
--- a/src/nrniv/nrndae.cpp
+++ b/src/nrniv/nrndae.cpp
@@ -4,6 +4,7 @@
 #include "nrndae_c.h"
 #include "nrnoc2iv.h"
 #include "treeset.h"
+#include "utils/enumerate.h"
 
 extern int secondorder;
 
@@ -24,16 +25,16 @@ void nrndae_deregister(NrnDAE* n) {
 
 int nrndae_extra_eqn_count() {
     int neqn = 0;
-    for (NrnDAEPtrListIterator m = nrndae_list.begin(); m != nrndae_list.end(); m++) {
-        neqn += (*m)->extra_eqn_count();
+    for (NrnDAE* item: nrndae_list) {
+        neqn += item->extra_eqn_count();
     }
     return neqn;
 }
 
 void nrndae_update(NrnThread* _nt) {
     update_sp13_rhs_based_on_actual_rhs(_nt);
-    for (NrnDAEPtrListIterator m = nrndae_list.begin(); m != nrndae_list.end(); m++) {
-        (*m)->update();
+    for (NrnDAE* item: nrndae_list) {
+        item->update();
     }
     update_actual_rhs_based_on_sp13_rhs(_nt);
 }
@@ -45,9 +46,9 @@ void nrndae_alloc() {
     if (_nt->_ecell_memb_list) {
         neqn += _nt->_ecell_memb_list->nodecount * nlayer;
     }
-    for (NrnDAEPtrListIterator m = nrndae_list.begin(); m != nrndae_list.end(); m++) {
-        (*m)->alloc(neqn + 1);
-        neqn += (*m)->extra_eqn_count();
+    for (NrnDAE* item: nrndae_list) {
+        item->alloc(neqn + 1);
+        neqn += item->extra_eqn_count();
     }
 }
 
@@ -62,8 +63,8 @@ void nrndae_init() {
         (secondorder > 0 || ((cvode_active_ > 0) && (nrn_use_daspk_ == 0)))) {
         hoc_execerror("NrnDAEs only work with secondorder==0 or daspk", 0);
     }
-    for (NrnDAEPtrListIterator m = nrndae_list.begin(); m != nrndae_list.end(); m++) {
-        (*m)->init();
+    for (NrnDAE* item: nrndae_list) {
+        item->init();
     }
     for (int it = 0; it < nrn_nthread; ++it) {
         auto* const nt = std::next(nrn_threads, it);
@@ -75,31 +76,31 @@ void nrndae_init() {
 void nrndae_rhs(NrnThread* _nt) {
     update_sp13_mat_based_on_actual_d(_nt);
     update_sp13_rhs_based_on_actual_rhs(_nt);
-    for (NrnDAEPtrListIterator m = nrndae_list.begin(); m != nrndae_list.end(); m++) {
-        (*m)->rhs();
+    for (NrnDAE* item: nrndae_list) {
+        item->rhs();
     }
     update_actual_d_based_on_sp13_mat(_nt);
     update_actual_rhs_based_on_sp13_rhs(_nt);
 }
 
 void nrndae_lhs() {
-    for (NrnDAEPtrListIterator m = nrndae_list.begin(); m != nrndae_list.end(); m++) {
-        (*m)->lhs();
+    for (NrnDAE* item: nrndae_list) {
+        item->lhs();
     }
 }
 
 void nrndae_dkmap(std::vector<neuron::container::data_handle<double>>& pv,
                   std::vector<neuron::container::data_handle<double>>& pvdot) {
-    for (NrnDAEPtrListIterator m = nrndae_list.begin(); m != nrndae_list.end(); m++) {
-        (*m)->dkmap(pv, pvdot);
+    for (NrnDAE* item: nrndae_list) {
+        item->dkmap(pv, pvdot);
     }
 }
 
 void nrndae_dkres(double* y, double* yprime, double* delta) {
     // c*y' = f(y) so
     // delta = c*y' - f(y)
-    for (NrnDAEPtrListIterator m = nrndae_list.begin(); m != nrndae_list.end(); m++) {
-        (*m)->dkres(y, yprime, delta);
+    for (NrnDAE* item: nrndae_list) {
+        item->dkres(y, yprime, delta);
     }
 }
 
diff --git a/src/nrniv/prcellstate.cpp b/src/nrniv/prcellstate.cpp
index eb7b5e4822..46e4a0952b 100644
--- a/src/nrniv/prcellstate.cpp
+++ b/src/nrniv/prcellstate.cpp
@@ -6,6 +6,7 @@
 #include <map>
 #include "OS/list.h"
 #include "neuron.h"
+#include "utils/enumerate.h"
 
 #define precision 15
 
@@ -47,11 +48,7 @@ static void pr_netcon(NrnThread& nt, FILE* f, const std::map<void*, int>& pnt2in
 
     // List of NetCon for each of the NET_RECEIVE point process instances
     // ... all NetCon list in the hoc NetCon cTemplate
-    std::vector<NetCon*>** nclist = new std::vector<NetCon*>*[pnt2index.size()];
-    for (size_t i = 0; i < pnt2index.size(); ++i) {
-        nclist[i] = new std::vector<NetCon*>();
-        nclist[i]->reserve(1);
-    }
+    std::vector<std::vector<NetCon*>> nclist(pnt2index.size());
     int nc_cnt = 0;
     Symbol* ncsym = hoc_lookup("NetCon");
     hoc_Item* q;
@@ -61,17 +58,15 @@ static void pr_netcon(NrnThread& nt, FILE* f, const std::map<void*, int>& pnt2in
         Point_process* pp = nc->target_;
         const auto& it = pnt2index.find(pp);
         if (it != pnt2index.end()) {
-            nclist[it->second]->push_back(nc);
+            nclist[it->second].push_back(nc);
             ++nc_cnt;
         }
     }
     fprintf(f, "netcons %d\n", nc_cnt);
     fprintf(f, " pntindex srcgid active delay weights\n");
-    for (size_t i = 0; i < pnt2index.size(); ++i) {
-        for (int j = 0; j < nclist[i]->size(); ++j) {
-            NetCon* nc = nclist[i]->at(j);
-            int srcgid = -3;
-            srcgid = (nc->src_) ? nc->src_->gid_ : -3;
+    for (const auto&& [i, ncl]: enumerate(nclist)) {
+        for (const auto& nc: ncl) {
+            int srcgid = (nc->src_) ? nc->src_->gid_ : -3;
             if (srcgid < 0 && nc->src_ && nc->src_->osrc_) {
                 const char* name = nc->src_->osrc_->ctemplate->sym->name;
                 fprintf(f, "%zd %s %d %.*g", i, name, nc->active_ ? 1 : 0, precision, nc->delay_);
@@ -87,11 +82,6 @@ static void pr_netcon(NrnThread& nt, FILE* f, const std::map<void*, int>& pnt2in
             fprintf(f, "\n");
         }
     }
-    // cleanup
-    for (size_t i = 0; i < pnt2index.size(); ++i) {
-        delete nclist[i];
-    }
-    delete[] nclist;
 }
 
 static void pr_realcell(PreSyn& ps, NrnThread& nt, FILE* f) {
diff --git a/src/nrniv/pysecname2sec.cpp b/src/nrniv/pysecname2sec.cpp
index a47be3e2f4..2f0d4a0201 100644
--- a/src/nrniv/pysecname2sec.cpp
+++ b/src/nrniv/pysecname2sec.cpp
@@ -8,6 +8,7 @@
 #include <nrnoc2iv.h>
 
 #include "nrnsymdiritem.h"
+#include "utils/enumerate.h"
 
 #include <string>
 #include <map>
@@ -249,10 +250,9 @@ void nrn_symdir_load_pysec(std::vector<SymbolItem*>& sl, void* v) {
     activate();
     if (!v) {
         // top level items are any of the four types
-        for (Name2CellorSec::iterator it = n2cs.begin(); it != n2cs.end(); ++it) {
-            CellorSec& cs = it->second;
+        for (auto&& [symbol, cs]: n2cs) {
             if (cs.first != NONETYPE && cs.first != OVERLOADCOUNT) {
-                SymbolItem* si = new SymbolItem(it->first.c_str(), 0);
+                SymbolItem* si = new SymbolItem(symbol.c_str(), 0);
                 si->pysec_type_ = cs.first == CELLTYPE ? PYSECOBJ : PYSECNAME;
                 si->pysec_ = (Section*) cs.second;
                 sl.push_back(si);
@@ -260,11 +260,9 @@ void nrn_symdir_load_pysec(std::vector<SymbolItem*>& sl, void* v) {
         }
     } else {
         // in cell items are either OVERLOADCOUNT or SECTYPE
-        Name2CellorSec* n2s = (Name2CellorSec*) v;
-        for (Name2CellorSec::iterator it = n2s->begin(); it != n2s->end(); ++it) {
-            CellorSec& cs = it->second;
+        for (auto&& [symbol, cs]: *static_cast<Name2CellorSec*>(v)) {
             if (cs.first == SECTYPE) {
-                SymbolItem* si = new SymbolItem(it->first.c_str(), 0);
+                auto* si = new SymbolItem(symbol.c_str(), 0);
                 si->pysec_type_ = PYSECNAME;
                 si->pysec_ = (Section*) cs.second;
                 sl.push_back(si);
diff --git a/src/nrniv/savstate.cpp b/src/nrniv/savstate.cpp
index 0e2df0f2dc..eeb6ca81bd 100644
--- a/src/nrniv/savstate.cpp
+++ b/src/nrniv/savstate.cpp
@@ -13,6 +13,7 @@
 #include "tqueue.h"
 #include "netcon.h"
 #include "vrecitem.h"
+#include "utils/enumerate.h"
 
 typedef void (*ReceiveFunc)(Point_process*, double*, double);
 
@@ -606,11 +607,10 @@ void SaveState::save() {
         assert(t == nt->_t);
     }
     t_ = t;
-    int inode;
     for (int isec = 0; isec < nsec_; ++isec) {
         SecState& ss = ss_[isec];
         Section* sec = ss.sec;
-        for (inode = 0; inode < ss.nnode; ++inode) {
+        for (int inode = 0; inode < ss.nnode; ++inode) {
             NodeState& ns = ss.ns[inode];
             Node* nd = sec->pnode[inode];
             savenode(ns, nd);
@@ -629,8 +629,8 @@ void SaveState::save() {
     if (nprs_) {
         std::vector<PlayRecord*>* prl = net_cvode_instance_prl();
         assert(nprs_ == prl->size());
-        for (std::size_t i = 0; i < nprs_; ++i) {
-            prs_[i] = (*prl)[i]->savestate_save();
+        for (auto&& [i, e]: enumerate(*prl)) {
+            prs_[i] = e->savestate_save();
         }
     }
     savenet();
@@ -645,8 +645,7 @@ void SaveState::save() {
 void SaveState::savenode(NodeState& ns, Node* nd) {
     ns.v = NODEV(nd);
     int istate = 0;
-    Prop* p;
-    for (p = nd->prop; p; p = p->next) {
+    for (Prop* p = nd->prop; p; p = p->next) {
         if (ssi[p->_type].size == 0) {
             continue;
         }
@@ -687,11 +686,10 @@ void SaveState::restore(int type) {
     FOR_THREADS(nt) {
         nt->_t = t_;
     }
-    int inode;
     for (int isec = 0; isec < nsec_; ++isec) {
         SecState& ss = ss_[isec];
         Section* sec = ss.sec;
-        for (inode = 0; inode < ss.nnode; ++inode) {
+        for (int inode = 0; inode < ss.nnode; ++inode) {
             NodeState& ns = ss.ns[inode];
             Node* nd = sec->pnode[inode];
             restorenode(ns, nd);
@@ -714,8 +712,7 @@ void SaveState::restore(int type) {
     // during a local step prl is augmented with GLineRecord
     // assert(nprs_ == prl->count());
     assert(nprs_ <= prl->size());
-    int i;
-    for (i = 0; i < nprs_; ++i) {
+    for (int i = 0; i < nprs_; ++i) {
         prs_[i]->savestate_restore();
     }
     restorenet();
@@ -730,8 +727,7 @@ void SaveState::restore(int type) {
 void SaveState::restorenode(NodeState& ns, Node* nd) {
     nd->v() = ns.v;
     int istate = 0;
-    Prop* p;
-    for (p = nd->prop; p; p = p->next) {
+    for (Prop* p = nd->prop; p; p = p->next) {
         if (ssi[p->_type].size == 0) {
             continue;
         }
@@ -933,27 +929,21 @@ void SaveState::write(OcFile* ocf, bool close) {
 }
 
 void SaveState::savenet() {
-    int i, j, n;
-    double* w;
     hoc_Item* q;
-    Object* ob;
-    NetCon* d;
-    PreSyn* ps;
-    i = 0;
+    int i = 0;
     ITERATE(q, nct->olist) {
-        ob = OBJ(q);
-        d = (NetCon*) ob->u.this_pointer;
-        n = ncs_[i].nstate;
-        w = ncs_[i].state;
-        for (j = 0; j < n; ++j) {
+        Object* ob = OBJ(q);
+        const NetCon* d = (NetCon*) ob->u.this_pointer;
+        int n = ncs_[i].nstate;
+        double* w = ncs_[i].state;
+        for (int j = 0; j < n; ++j) {
             w[j] = d->weight_[j];
         }
         ++i;
     }
-    i = 0;
-    if (net_cvode_instance_psl())
+    if (int i = 0; net_cvode_instance_psl()) {
         ITERATE(q, net_cvode_instance_psl()) {
-            ps = (PreSyn*) VOIDITM(q);
+            auto* ps = static_cast<PreSyn*>(VOIDITM(q));
             ps->hi_index_ = i;
             pss_[i].flag = ps->flag_;
             pss_[i].valthresh = ps->valthresh_;
@@ -961,6 +951,7 @@ void SaveState::savenet() {
             pss_[i].told = ps->told_;
             ++i;
         }
+    }
     alloc_tq();
     tqcnt_ = 0;
     NrnThread* nt;
@@ -984,29 +975,23 @@ void SaveState::tqsave(const TQItem* q, int) {
 }
 
 void SaveState::restorenet() {
-    int i, j, n;
-    double* w;
-    hoc_Item* q;
-    Object* ob;
-    NetCon* d;
-    PreSyn* ps;
     // NetCon's
-    i = 0;
+    int i = 0;
+    hoc_Item* q;
     ITERATE(q, nct->olist) {
-        ob = OBJ(q);
-        d = (NetCon*) ob->u.this_pointer;
-        n = ncs_[i].nstate;
-        w = ncs_[i].state;
-        for (j = 0; j < n; ++j) {
+        Object* ob = OBJ(q);
+        NetCon* d = (NetCon*) ob->u.this_pointer;
+        int n = ncs_[i].nstate;
+        const double* w = ncs_[i].state;
+        for (int j = 0; j < n; ++j) {
             d->weight_[j] = w[j];
         }
         ++i;
     }
     // PreSyn's
-    i = 0;
-    if (net_cvode_instance_psl())
+    if (int i = 0; net_cvode_instance_psl())
         ITERATE(q, net_cvode_instance_psl()) {
-            ps = (PreSyn*) VOIDITM(q);
+            auto* ps = static_cast<PreSyn*>(VOIDITM(q));
             ps->hi_index_ = i;
             ps->flag_ = pss_[i].flag;
             ps->valthresh_ = pss_[i].valthresh;
@@ -1019,8 +1004,7 @@ void SaveState::restorenet() {
     // clear it
     clear_event_queue();
     // restore it
-    n = tqs_->nstate;
-    for (i = 0; i < n; ++i) {
+    for (int i = 0; i < tqs_->nstate; ++i) {
         tqs_->items[i]->savestate_restore(tqs_->tdeliver[i], net_cvode_instance);
     }
 }
@@ -1033,8 +1017,7 @@ void SaveState::readnet(FILE* f) {
     if (nncs_ != 0) {
         ncs_ = new NetConState[nncs_];
     }
-    int i, n, type;
-    for (i = 0; i < nncs_; ++i) {
+    for (int i = 0; i < nncs_; ++i) {
         ASSERTfgets(buf, 200, f);
         sscanf(buf, "%d %d\n", &ncs_[i].object_index, &ncs_[i].nstate);
         if (ncs_[i].nstate) {
@@ -1049,7 +1032,7 @@ void SaveState::readnet(FILE* f) {
         pss_ = new PreSynState[npss_];
         ASSERTfread(pss_, sizeof(PreSynState), npss_, f);
         PreSyn* ps;
-        i = 0;
+        int i = 0;
         hoc_Item* q;
         if (net_cvode_instance_psl())
             ITERATE(q, net_cvode_instance_psl()) {
@@ -1060,6 +1043,7 @@ void SaveState::readnet(FILE* f) {
         assert(npss_ == i);
     }
 
+    int n = 0;
     ASSERTfgets(buf, 200, f);
     sscanf(buf, "%d\n", &n);
     tqs_->nstate = n;
@@ -1067,9 +1051,10 @@ void SaveState::readnet(FILE* f) {
         tqs_->items = new DiscreteEvent*[n];
         tqs_->tdeliver = new double[n];
         ASSERTfread(tqs_->tdeliver, sizeof(double), n, f);
-        for (i = 0; i < n; ++i) {
+        for (int i = 0; i < n; ++i) {
             DiscreteEvent* de = NULL;
             ASSERTfgets(buf, 200, f);
+            int type = 0;
             sscanf(buf, "%d\n", &type);
             switch (type) {
             case DiscreteEventType:
@@ -1104,8 +1089,7 @@ void SaveState::readnet(FILE* f) {
 
 void SaveState::writenet(FILE* f) {
     fprintf(f, "%d\n", nncs_);
-    int i, n;
-    for (i = 0; i < nncs_; ++i) {
+    for (int i = 0; i < nncs_; ++i) {
         fprintf(f, "%d %d\n", ncs_[i].object_index, ncs_[i].nstate);
         if (ncs_[i].nstate) {
             ASSERTfwrite(ncs_[i].state, sizeof(double), ncs_[i].nstate, f);
@@ -1115,11 +1099,11 @@ void SaveState::writenet(FILE* f) {
     if (npss_) {
         ASSERTfwrite(pss_, sizeof(PreSynState), npss_, f);
     }
-    n = tqs_->nstate;
+    int n = tqs_->nstate;
     fprintf(f, "%d\n", n);
     if (n) {
         ASSERTfwrite(tqs_->tdeliver, sizeof(double), n, f);
-        for (i = 0; i < n; ++i) {
+        for (int i = 0; i < n; ++i) {
             tqs_->items[i]->savestate_write(f);
         }
     }
@@ -1135,14 +1119,11 @@ bool SaveState::checknet(bool warn) {
         }
         return false;
     }
-    int i;
     hoc_Item* q;
-    Object* ob;
-    NetCon* d;
-    i = 0;
+    int i = 0;
     ITERATE(q, nct->olist) {
-        ob = OBJ(q);
-        d = (NetCon*) ob->u.this_pointer;
+        Object* ob = OBJ(q);
+        const auto* d = static_cast<NetCon*>(ob->u.this_pointer);
         if (ob->index != ncs_[i].object_index) {
             if (warn) {
                 fprintf(stderr,
@@ -1187,14 +1168,11 @@ void SaveState::allocnet() {
     if (nncs_ != 0) {
         ncs_ = new NetConState[nncs_];
     }
-    int i, n;
     hoc_Item* q;
-    Object* ob;
-    NetCon* d;
-    i = 0;
+    int i = 0;
     ITERATE(q, nct->olist) {
-        ob = OBJ(q);
-        d = (NetCon*) ob->u.this_pointer;
+        Object* ob = OBJ(q);
+        const auto* d = static_cast<NetCon*>(ob->u.this_pointer);
         ncs_[i].object_index = ob->index;
         ncs_[i].nstate = d->cnt_;
         if (d->cnt_) {
@@ -1202,11 +1180,10 @@ void SaveState::allocnet() {
         }
         ++i;
     }
-    PreSyn* ps;
     npss_ = 0;
     if (net_cvode_instance_psl())
         ITERATE(q, net_cvode_instance_psl()) {
-            ps = (PreSyn*) VOIDITM(q);
+            auto* ps = static_cast<PreSyn*>(VOIDITM(q));
             ps->hi_index_ = npss_;
             ++npss_;
         }
@@ -1218,9 +1195,8 @@ void SaveState::allocnet() {
 // The event TQueue is highly volatile so it needs to be freed and allocated
 // on every save and fread
 void SaveState::free_tq() {
-    int i;
     if (tqs_->nstate) {
-        for (i = 0; i < tqs_->nstate; ++i) {
+        for (int i = 0; i < tqs_->nstate; ++i) {
             delete tqs_->items[i];
         }
         tqs_->nstate = 0;
@@ -1229,7 +1205,6 @@ void SaveState::free_tq() {
     }
 }
 void SaveState::alloc_tq() {
-    int n;
     free_tq();
     tqcnt_ = 0;
     NrnThread* nt;
@@ -1239,7 +1214,7 @@ void SaveState::alloc_tq() {
         callback_mode = 0;
         tq->forall_callback(tqcallback);
     }
-    n = tqcnt_;
+    int n = tqcnt_;
     tqs_->nstate = n;
     if (n) {
         tqs_->items = new DiscreteEvent*[n];
diff --git a/src/oc/code.cpp b/src/oc/code.cpp
index 6aa92e2e4c..d00a24bf7d 100644
--- a/src/oc/code.cpp
+++ b/src/oc/code.cpp
@@ -17,6 +17,7 @@
 #include <stdlib.h>
 #include <nrnmpi.h>
 #include "nrnfilewrap.h"
+#include "utils/enumerate.h"
 
 
 #include "options.h"
@@ -398,15 +399,13 @@ static void frame_objauto_recover_on_err(Frame* ff) { /* only on error */
 static void stack_obtmp_recover_on_err(int tcnt) {
     if (tobj_count > tcnt) {
         // unref tmpobjects from the top of the stack until we have the right number left
-        for (auto stkp = stack.rbegin(); stkp != stack.rend(); ++stkp) {
-            // What is the index in `stack` of `stkp`?
-            auto const index = stack.size() - 1 - std::distance(stack.rbegin(), stkp);
-            if (stack_entry_is_tmpobject(*stkp)) {
-                hoc_stkobj_unref(std::get<Object*>(*stkp), index);
+        for (const auto&& [index, stkp]: renumerate(stack)) {
+            if (stack_entry_is_tmpobject(stkp)) {
+                hoc_stkobj_unref(std::get<Object*>(stkp), index);
                 if (tobj_count == tcnt) {
                     return;
                 }
-            } else if (std::holds_alternative<std::nullptr_t>(*stkp)) {
+            } else if (std::holds_alternative<std::nullptr_t>(stkp)) {
                 printf("OBJECTTMP at stack index %ld already unreffed\n", index);
             }
         }
@@ -439,7 +438,7 @@ void hoc_prstack() {
     std::size_t i{};
     std::ostringstream oss;
     oss << "interpreter stack: " << stack.size() << '\n';
-    for (auto stkp = stack.rbegin(); stkp != stack.rend(); ++stkp, ++i) {
+    for (auto&& stkp: reverse(stack)) {
         if (i > 10) {
             oss << " ...\n";
             break;
@@ -454,7 +453,8 @@ void hoc_prstack() {
                 }
                 oss << ' ' << cxx_demangle(typeid(decltype(value)).name()) << '\n';
             },
-            *stkp);
+            stkp);
+        ++i;
     }
     Printf(oss.str().c_str());
 }
diff --git a/src/utils/enumerate.h b/src/utils/enumerate.h
new file mode 100644
index 0000000000..1a0b38f2ad
--- /dev/null
+++ b/src/utils/enumerate.h
@@ -0,0 +1,145 @@
+#pragma once
+
+#include <algorithm>
+#include <functional>
+#include <tuple>
+
+template <typename T,
+          typename F,
+          typename value_type = typename T::value_type,
+          typename TIter = decltype(std::begin(std::declval<T>())),
+          typename = decltype(std::end(std::declval<T>()))>
+void apply_to_first(T&& iterable, value_type&& value, F&& f) {
+    auto it = std::find(std::begin(std::forward<T>(iterable)),
+                        std::end(std::forward<T>(iterable)),
+                        std::forward<value_type>(value));
+    if (it != std::end(std::forward<T>(iterable))) {
+        f(it);
+    }
+}
+
+template <typename T, typename value_type = typename T::value_type>
+void erase_first(T&& iterable, value_type&& value) {
+    apply_to_first(std::forward<T>(iterable),
+                   std::forward<value_type>(value),
+                   [&iterable](const auto& it) { iterable.erase(it); });
+}
+
+template <typename T,
+          typename TIter = decltype(std::begin(std::declval<T>())),
+          typename = decltype(std::begin(std::declval<T>())),
+          typename = decltype(std::end(std::declval<T>()))>
+constexpr auto range(T&& iterable) {
+    struct iterator {
+        std::size_t i;
+        TIter iter;
+        bool operator!=(const iterator& other) const {
+            return iter != other.iter;
+        }
+        void operator++() {
+            ++i;
+            ++iter;
+        }
+        auto operator*() const {
+            return i;
+        }
+    };
+    struct iterable_wrapper {
+        T iterable;
+        auto begin() {
+            return iterator{0, std::begin(iterable)};
+        }
+        auto end() {
+            return iterator{0, std::end(iterable)};
+        }
+    };
+    return iterable_wrapper{std::forward<T>(iterable)};
+}
+
+template <typename T,
+          typename TIter = decltype(std::rbegin(std::declval<T>())),
+          typename = decltype(std::rend(std::declval<T>()))>
+constexpr auto reverse(T&& iterable) {
+    struct iterator {
+        TIter iter;
+        bool operator!=(const iterator& other) const {
+            return iter != other.iter;
+        }
+        void operator++() {
+            ++iter;
+        }
+        auto&& operator*() const {
+            return *iter;
+        }
+    };
+    struct iterable_wrapper {
+        T iterable;
+        auto begin() {
+            return iterator{std::rbegin(iterable)};
+        }
+        auto end() {
+            return iterator{std::rend(iterable)};
+        }
+    };
+    return iterable_wrapper{std::forward<T>(iterable)};
+}
+
+template <typename T,
+          typename TIter = decltype(std::begin(std::declval<T>())),
+          typename = decltype(std::end(std::declval<T>()))>
+constexpr auto enumerate(T&& iterable) {
+    struct iterator {
+        std::size_t i;
+        TIter iter;
+        bool operator!=(const iterator& other) const {
+            return iter != other.iter;
+        }
+        void operator++() {
+            ++i;
+            ++iter;
+        }
+        auto operator*() const {
+            return std::tie(i, *iter);
+        }
+    };
+    struct iterable_wrapper {
+        T iterable;
+        auto begin() {
+            return iterator{0, std::begin(iterable)};
+        }
+        auto end() {
+            return iterator{0, std::end(iterable)};
+        }
+    };
+    return iterable_wrapper{std::forward<T>(iterable)};
+}
+
+template <typename T,
+          typename TIter = decltype(std::rbegin(std::declval<T>())),
+          typename = decltype(std::rend(std::declval<T>()))>
+constexpr auto renumerate(T&& iterable) {
+    struct iterator {
+        std::size_t i;
+        TIter iter;
+        bool operator!=(const iterator& other) const {
+            return iter != other.iter;
+        }
+        void operator++() {
+            --i;
+            ++iter;
+        }
+        auto operator*() const {
+            return std::tie(i, *iter);
+        }
+    };
+    struct iterable_wrapper {
+        T iterable;
+        auto begin() {
+            return iterator{std::size(iterable) - 1, std::rbegin(iterable)};
+        }
+        auto end() {
+            return iterator{std::size(iterable) - 1, std::rend(iterable)};
+        }
+    };
+    return iterable_wrapper{std::forward<T>(iterable)};
+}
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 0fb11833a2..729a95d771 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -22,6 +22,7 @@ add_executable(
   unit_tests/container/generic_data_handle.cpp
   unit_tests/container/mechanism.cpp
   unit_tests/container/node.cpp
+  unit_tests/utils/enumerate.cpp
   unit_tests/oc/hoc_interpreter.cpp)
 set(catch2_targets testneuron)
 if(NRN_ENABLE_THREADS)
diff --git a/test/unit_tests/utils/enumerate.cpp b/test/unit_tests/utils/enumerate.cpp
new file mode 100644
index 0000000000..2ef2dd00bf
--- /dev/null
+++ b/test/unit_tests/utils/enumerate.cpp
@@ -0,0 +1,85 @@
+#include <vector>
+
+#include "utils/enumerate.h"
+
+#include <catch2/catch.hpp>
+
+
+TEST_CASE("apply_to_first", "[Neuron]") {
+    std::vector<double> x{1.0, 2.0, 2.0, 3.0};
+
+    apply_to_first(x, 2.0, [](auto it) { *it = 5.0; });
+    REQUIRE(x == std::vector<double>({1.0, 5.0, 2.0, 3.0}));
+}
+
+TEST_CASE("erase_first", "[Neuron]") {
+    std::vector<double> x{1.0, 2.0, 2.0, 3.0};
+
+    erase_first(x, 2.0);
+    REQUIRE(x == std::vector<double>({1.0, 2.0, 3.0}));
+}
+
+TEST_CASE("reverse", "[Neuron]") {
+    std::vector<double> x{1.0, 2.0, 3.0};
+
+    for (auto& i: reverse(x)) {
+        i *= -1.0;
+    }
+    REQUIRE(x == std::vector<double>({-1.0, -2.0, -3.0}));
+}
+
+TEST_CASE("reverse; no-copy", "[Neuron]") {
+    std::vector<double> x{1.0, 2.0, 3.0};
+
+    auto reverse_iterable = reverse(x);
+
+    for (auto& xx: x) {
+        xx *= -1.0;
+    }
+
+    for (const auto& xx: reverse_iterable) {
+        REQUIRE(xx < 0.0);
+    }
+}
+
+TEST_CASE("range", "[Neuron]") {
+    std::vector<double> x{1.0, 2.0, 3.0};
+
+    std::vector<std::size_t> v{};
+    for (std::size_t i: range(x)) {
+        v.push_back(i);
+    }
+    REQUIRE(v == std::vector<size_t>{0, 1, 2});
+}
+
+TEST_CASE("enumerate", "[Neuron]") {
+    std::vector<double> x{1.0, 2.0, 3.0};
+
+    int j = 0;
+    for (auto&& [i, elem]: enumerate(x)) {
+        if (i == 0)
+            REQUIRE(elem == 1.0);
+        if (i == 1)
+            REQUIRE(elem == 2.0);
+        if (i == 2)
+            REQUIRE(elem == 3.0);
+        REQUIRE(i == j);
+        ++j;
+    }
+}
+
+TEST_CASE("renumerate", "[Neuron]") {
+    std::vector<double> x{1.0, 2.0, 3.0};
+
+    int j = x.size() - 1;
+    for (auto&& [i, elem]: renumerate(x)) {
+        if (i == 0)
+            REQUIRE(elem == 1.0);
+        if (i == 1)
+            REQUIRE(elem == 2.0);
+        if (i == 2)
+            REQUIRE(elem == 3.0);
+        REQUIRE(i == j);
+        --j;
+    }
+}