diff --git a/.circleci/config.yml b/.circleci/config.yml
index 263a46c424..c1d872ad41 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -41,6 +41,7 @@ jobs:
       - run:
           name: Test manylinux AArch64 wheel
           command: |
+            sudo apt purge -y needrestart
             sudo apt update
             echo "deb http://ppa.launchpad.net/deadsnakes/ppa/ubuntu jammy main" | sudo tee /etc/apt/sources.list.d/deadsnakes-ppa.list && sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-keys F23C5A6CF475977595C89F51BA6932366A755776 && sudo apt update
             sudo apt update
diff --git a/src/nrncvode/cvodeobj.cpp b/src/nrncvode/cvodeobj.cpp
index cdb05b511e..89b2ca3e58 100644
--- a/src/nrncvode/cvodeobj.cpp
+++ b/src/nrncvode/cvodeobj.cpp
@@ -956,14 +956,12 @@ void Cvode::maxstate(bool b, NrnThread* nt) {
 }
 
 void Cvode::maxstate(double* pd) {
-    int i;
-    NrnThread* nt;
     if (maxstate_) {
-        FOR_THREADS(nt) {
+        for (NrnThread* nt: for_threads(nrn_threads, nrn_nthread)) {
             double* m = n_vector_data(maxstate_, nt->id);
             int n = ctd_[nt->id].nvsize_;
             int o = ctd_[nt->id].nvoffset_;
-            for (i = 0; i < n; ++i) {
+            for (int i = 0; i < n; ++i) {
                 pd[i + o] = m[i];
             }
         }
@@ -971,14 +969,12 @@ void Cvode::maxstate(double* pd) {
 }
 
 void Cvode::maxacor(double* pd) {
-    int i;
-    NrnThread* nt;
     if (maxacor_) {
-        FOR_THREADS(nt) {
+        for (const NrnThread* nt: for_threads(nrn_threads, nrn_nthread)) {
             double* m = n_vector_data(maxacor_, nt->id);
             int n = ctd_[nt->id].nvsize_;
             int o = ctd_[nt->id].nvoffset_;
-            for (i = 0; i < n; ++i) {
+            for (int i = 0; i < n; ++i) {
                 pd[i + o] = m[i];
             }
         }
@@ -1221,13 +1217,12 @@ int Cvode::init(double tout) {
 }
 
 int Cvode::interpolate(double tout) {
-    NrnThread* _nt;
     if (neq_ == 0) {
         t_ = tout;
         if (nth_) {
             nth_->_t = t_;
         } else {
-            FOR_THREADS(_nt) {
+            for (NrnThread* _nt: for_threads(nrn_threads, nrn_nthread)) {
                 _nt->_t = t_;
             }
         }
@@ -1241,7 +1236,7 @@ int Cvode::interpolate(double tout) {
         if (nth_) {  // lvardt
             nth_->_t = tout;
         } else {
-            FOR_THREADS(_nt) {
+            for (NrnThread* _nt: for_threads(nrn_threads, nrn_nthread)) {
                 _nt->_t = tout;  // but leave t_ at the initialization point.
             }
         }
diff --git a/src/nrncvode/netcvode.cpp b/src/nrncvode/netcvode.cpp
index 7df7471b4d..0110299653 100644
--- a/src/nrncvode/netcvode.cpp
+++ b/src/nrncvode/netcvode.cpp
@@ -1460,7 +1460,6 @@ bool NetCvode::init_global() {
     structure_change_cnt_ = structure_change_cnt;
     matrix_change_cnt_ = -1;
     playrec_change_cnt_ = 0;
-    NrnThread* _nt;
     // We copy Memb_list* into cml->ml below. At the moment this CVode code
     // generates its own complicated set of Memb_list* that operate in
     // list-of-handles mode instead of referring to contiguous sets of values.
@@ -1486,7 +1485,7 @@ bool NetCvode::init_global() {
         del_cv_memb_list();
         Cvode& cv = *gcv_;
         distribute_dinfo(nullptr, 0);
-        FOR_THREADS(_nt) {
+        for (NrnThread* _nt: for_threads(nrn_threads, nrn_nthread)) {
             CvodeThreadData& z = cv.ctd_[_nt->id];
             z.rootnodecount_ = _nt->ncell;
             z.v_node_count_ = _nt->end;
@@ -3963,7 +3962,6 @@ void NetCvode::re_init(double t) {
 }
 
 void NetCvode::fornetcon_prepare() {
-    NrnThread* nt;
     NrnThreadMembList* tml;
     if (fornetcon_change_cnt_ == structure_change_cnt) {
         return;
@@ -3994,17 +3992,19 @@ void NetCvode::fornetcon_prepare() {
                 fnc->size = 0;
             }
         } else {
-            FOR_THREADS(nt) for (tml = nt->tml; tml; tml = tml->next) if (tml->index == type) {
-                Memb_list* m = tml->ml;
-                for (j = 0; j < m->nodecount; ++j) {
-                    void** v = &(m->pdata[j][index].literal_value<void*>());
-                    _nrn_free_fornetcon(v);
-                    ForNetConsInfo* fnc = new ForNetConsInfo;
-                    *v = fnc;
-                    fnc->argslist = 0;
-                    fnc->size = 0;
-                }
-            }
+            for (NrnThread* nt: for_threads(nrn_threads, nrn_nthread))
+                for (tml = nt->tml; tml; tml = tml->next)
+                    if (tml->index == type) {
+                        Memb_list* m = tml->ml;
+                        for (j = 0; j < m->nodecount; ++j) {
+                            void** v = &(m->pdata[j][index].literal_value<void*>());
+                            _nrn_free_fornetcon(v);
+                            ForNetConsInfo* fnc = new ForNetConsInfo;
+                            *v = fnc;
+                            fnc->argslist = 0;
+                            fnc->size = 0;
+                        }
+                    }
         }
     }
     // two loops over all netcons. one to count, one to fill in argslist
@@ -4037,18 +4037,19 @@ void NetCvode::fornetcon_prepare() {
                 }
             }
         } else {
-            FOR_THREADS(nt)
-            for (tml = nt->tml; tml; tml = tml->next)
-                if (tml->index == nrn_fornetcon_type_[i]) {
-                    Memb_list* m = tml->ml;
-                    for (j = 0; j < m->nodecount; ++j) {
-                        auto* fnc = static_cast<ForNetConsInfo*>(m->pdata[j][index].get<void*>());
-                        if (fnc->size > 0) {
-                            fnc->argslist = new double*[fnc->size];
-                            fnc->size = 0;
+            for (NrnThread* nt: for_threads(nrn_threads, nrn_nthread))
+                for (tml = nt->tml; tml; tml = tml->next)
+                    if (tml->index == nrn_fornetcon_type_[i]) {
+                        Memb_list* m = tml->ml;
+                        for (j = 0; j < m->nodecount; ++j) {
+                            auto* fnc = static_cast<ForNetConsInfo*>(
+                                m->pdata[j][index].get<void*>());
+                            if (fnc->size > 0) {
+                                fnc->argslist = new double*[fnc->size];
+                                fnc->size = 0;
+                            }
                         }
                     }
-                }
         }
     }
     // fill in argslist and count again
diff --git a/src/nrncvode/occvode.cpp b/src/nrncvode/occvode.cpp
index a2f6bd24eb..630e1d733e 100644
--- a/src/nrncvode/occvode.cpp
+++ b/src/nrncvode/occvode.cpp
@@ -99,7 +99,6 @@ bool Cvode::init_global() {
 void Cvode::init_eqn() {
     double vtol;
 
-    NrnThread* _nt;
     CvMembList* cml;
     int i, j, zneq, zneq_v, zneq_cap_v;
     // printf("Cvode::init_eqn\n");
@@ -125,7 +124,7 @@ void Cvode::init_eqn() {
         daspk_init_eqn();
         return;
     }
-    FOR_THREADS(_nt) {
+    for (NrnThread* _nt: for_threads(nrn_threads, nrn_nthread)) {
         // for lvardt, this body done only once and for ctd_[0]
         CvodeThreadData& z = ctd_[_nt->id];
         // how many ode's are there? First ones are non-zero capacitance
@@ -236,7 +235,7 @@ printf("%d Cvode::init_eqn id=%d neq_v_=%d #nonvint=%d #nonvint_extra=%d nvsize=
         z.no_cap_child_count_ = j;
 
         // use the sentinal values in NODERHS to construct a new no cap membrane list
-        new_no_cap_memb(z, _nt);
+        new_no_cap_memb(z, nullptr);
 
         // map the membrane mechanism ode state and dstate pointers
         int ieq = zneq_v;
@@ -270,7 +269,7 @@ printf("%d Cvode::init_eqn id=%d neq_v_=%d #nonvint=%d #nonvint_extra=%d nvsize=
     structure_change_ = false;
 }
 
-void Cvode::new_no_cap_memb(CvodeThreadData& z, NrnThread* _nt) {
+void Cvode::new_no_cap_memb(CvodeThreadData& z, NrnThread* /* thread */) {
     z.delete_memb_list(z.no_cap_memb_);
     z.no_cap_memb_ = nullptr;
     CvMembList* ncm{};
diff --git a/src/nrniv/nrncore_write/data/cell_group.cpp b/src/nrniv/nrncore_write/data/cell_group.cpp
index 8ee75d4191..85fb6f8607 100644
--- a/src/nrniv/nrncore_write/data/cell_group.cpp
+++ b/src/nrniv/nrncore_write/data/cell_group.cpp
@@ -552,9 +552,8 @@ void CellGroup::mk_tml_with_art(neuron::model_sorted_token const& cache_token, C
 size_t CellGroup::get_mla_rankbytes(CellGroup* cellgroups_) {
     size_t mla_rankbytes = 0;
     size_t nbytes;
-    NrnThread* nt;
     NrnThreadMembList* tml;
-    FOR_THREADS(nt) {
+    for (const NrnThread* nt: for_threads(nrn_threads, nrn_nthread)) {
         size_t threadbytes = 0;
         size_t npnt = 0;
         size_t nart = 0;
diff --git a/src/nrniv/savstate.cpp b/src/nrniv/savstate.cpp
index 190c935926..bc6e9f2b6a 100644
--- a/src/nrniv/savstate.cpp
+++ b/src/nrniv/savstate.cpp
@@ -594,11 +594,10 @@ void SaveState::ssfree() {
 }
 
 void SaveState::save() {
-    NrnThread* nt;
     if (!check(false)) {
         alloc();
     }
-    FOR_THREADS(nt) {
+    for (const NrnThread* nt: for_threads(nrn_threads, nrn_nthread)) {
         assert(t == nt->_t);
     }
     t_ = t;
@@ -673,12 +672,11 @@ void SaveState::saveacell(ACellState& ac, int type) {
 }
 
 void SaveState::restore(int type) {
-    NrnThread* nt;
     if (!check(true)) {
         hoc_execerror("SaveState:", "Stored state inconsistent with current neuron structure");
     }
     t = t_;
-    FOR_THREADS(nt) {
+    for (NrnThread* nt: for_threads(nrn_threads, nrn_nthread)) {
         nt->_t = t_;
     }
     for (int isec = 0; isec < nsec_; ++isec) {
@@ -948,8 +946,7 @@ void SaveState::savenet() {
     }
     alloc_tq();
     tqcnt_ = 0;
-    NrnThread* nt;
-    FOR_THREADS(nt) {
+    for (NrnThread* nt: for_threads(nrn_threads, nrn_nthread)) {
         TQueue* tq = net_cvode_instance_event_queue(nt);
         this_savestate = this;
         callback_mode = 1;
@@ -1195,8 +1192,7 @@ void SaveState::free_tq() {
 void SaveState::alloc_tq() {
     free_tq();
     tqcnt_ = 0;
-    NrnThread* nt;
-    FOR_THREADS(nt) {
+    for (NrnThread* nt: for_threads(nrn_threads, nrn_nthread)) {
         TQueue* tq = net_cvode_instance_event_queue(nt);
         this_savestate = this;
         callback_mode = 0;
diff --git a/src/nrnoc/extcelln.cpp b/src/nrnoc/extcelln.cpp
index ccfc4250a1..3ad4db7fa7 100644
--- a/src/nrnoc/extcelln.cpp
+++ b/src/nrnoc/extcelln.cpp
@@ -593,7 +593,7 @@ static void printnode(const char* s) {
 	Extnode* nde;	
 	double *pd;
 	NrnThread* _nt;
-	FOR_THREADS(_nt) for (in=0; in < _nt->end; ++in) {
+	for (NrnThread* _nt : for_threads(nrn_threads, nrn_nthread)) for (in=0; in < _nt->end; ++in) {
 		nd = _nt->_v_node[in];
 		if (nd->extnode) {
 			sec = nd->sec;
diff --git a/src/nrnoc/fadvance.cpp b/src/nrnoc/fadvance.cpp
index fea7150bae..c79348f709 100644
--- a/src/nrnoc/fadvance.cpp
+++ b/src/nrnoc/fadvance.cpp
@@ -804,7 +804,6 @@ void verify_structure(void) {
 
 void nrn_finitialize(int setv, double v) {
     int iord, i;
-    NrnThread* _nt;
     extern int _ninits;
     extern short* nrn_is_artificial_;
     ++_ninits;
@@ -834,14 +833,14 @@ void nrn_finitialize(int setv, double v) {
         nrn_deliver_events(nrn_threads + i); /* The play events at t=0 */
     }
     if (setv) {
-        FOR_THREADS(_nt) {
+        for (NrnThread* _nt: for_threads(nrn_threads, nrn_nthread)) {
             auto const vec_v = _nt->node_voltage_storage();
             std::fill_n(vec_v, _nt->end, v);
         }
     }
 #if 1 || NRNMPI
     if (nrnthread_vi_compute_)
-        FOR_THREADS(_nt) {
+        for (NrnThread* _nt: for_threads(nrn_threads, nrn_nthread)) {
             (*nrnthread_vi_compute_)(_nt);
         }
     {
@@ -850,7 +849,7 @@ void nrn_finitialize(int setv, double v) {
             (nrnmpi_v_transfer_)();
         }
         if (nrnthread_v_transfer_)
-            FOR_THREADS(_nt) {
+            for (NrnThread* _nt: for_threads(nrn_threads, nrn_nthread)) {
                 (*nrnthread_v_transfer_)(_nt);
             }
     }
diff --git a/src/nrnoc/multicore.cpp b/src/nrnoc/multicore.cpp
index 26d090248b..7827bce1b3 100644
--- a/src/nrnoc/multicore.cpp
+++ b/src/nrnoc/multicore.cpp
@@ -654,7 +654,6 @@ void nrn_thread_memblist_setup() {
 /* in passing, also set start and end indices. */
 
 void reorder_secorder() {
-    NrnThread* _nt;
     Section *sec, *ch;
     Node* nd;
     hoc_Item* qsec;
@@ -667,7 +666,7 @@ void reorder_secorder() {
         sec->order = -1;
     }
     order = 0;
-    FOR_THREADS(_nt) {
+    for (NrnThread* _nt: for_threads(nrn_threads, nrn_nthread)) {
         /* roots of this thread */
         sl = _nt->roots;
         inode = 0;
@@ -711,7 +710,7 @@ void reorder_secorder() {
         sec->order = -1;
     }
     order = 0;
-    FOR_THREADS(_nt) {
+    for (NrnThread* _nt: for_threads(nrn_threads, nrn_nthread)) {
         /* roots of this thread */
         sl = _nt->roots;
         inode = 0;
@@ -767,9 +766,10 @@ void reorder_secorder() {
       in either case, we can then point to v, d, rhs in proper
       node order
     */
-    FOR_THREADS(_nt) for (inode = 0; inode < _nt->end; ++inode) {
-        _nt->_v_node[inode]->_classical_parent = _nt->_v_parent[inode];
-    }
+    for (const NrnThread* _nt: for_threads(nrn_threads, nrn_nthread))
+        for (inode = 0; inode < _nt->end; ++inode) {
+            _nt->_v_node[inode]->_classical_parent = _nt->_v_parent[inode];
+        }
     if (nrn_multisplit_setup_) {
         /* classical order abandoned */
         (*nrn_multisplit_setup_)();
@@ -948,7 +948,6 @@ int nrn_user_partition() {
     hoc_List* sl;
     char buf[256];
     Section* sec;
-    NrnThread* nt;
     /* all one or all the other*/
     b = (nrn_threads[0].userpart != nullptr);
     for (it = 1; it < nrn_nthread; ++it) {
@@ -962,7 +961,7 @@ int nrn_user_partition() {
 
     /* discard partition if any section mentioned has been deleted. The
         model has changed */
-    FOR_THREADS(nt) {
+    for (NrnThread* nt: for_threads(nrn_threads, nrn_nthread)) {
         sl = nt->roots;
         ITERATE(qsec, sl) {
             sec = hocSEC(qsec);
@@ -983,7 +982,7 @@ int nrn_user_partition() {
     /* fill in ncell and verify consistency */
     n = 0;
     for (it = 0; it < nrn_nthread; ++it) {
-        nt = nrn_threads + it;
+        NrnThread* nt = nrn_threads + it;
         sl = nt->roots;
         nt->ncell = 0;
         ITERATE(qsec, sl) {
diff --git a/src/nrnoc/multicore.h b/src/nrnoc/multicore.h
index fb115b195c..dec977b9be 100644
--- a/src/nrnoc/multicore.h
+++ b/src/nrnoc/multicore.h
@@ -128,7 +128,43 @@ void reorder_secorder();
 void nrn_thread_memblist_setup();
 std::size_t nof_worker_threads();
 
-#define FOR_THREADS(nt) for (nt = nrn_threads; nt < nrn_threads + nrn_nthread; ++nt)
+
+// helper function for iterating over ``NrnThread``s
+inline auto for_threads(NrnThread* threads, int num_threads) {
+    struct iterator {
+        NrnThread* current;
+
+        NrnThread* operator*() const {
+            return current;
+        }
+        iterator& operator++() {
+            ++current;
+            return *this;
+        }
+        bool operator!=(const iterator& other) const {
+            return current != other.current;
+        }
+    };
+
+    struct iterable_wrapper {
+        NrnThread* base_;
+        int count_;
+
+        iterable_wrapper(NrnThread* base, int count)
+            : base_(base)
+            , count_(count) {}
+
+        iterator begin() const {
+            return iterator{base_};
+        }
+        iterator end() const {
+            return iterator{base_ + count_};
+        }
+    };
+
+    return iterable_wrapper(threads, num_threads);
+}
+
 
 // olupton 2022-01-31: could add a _NrnThread typedef here for .mod file
 //                     backwards compatibility if needed.
diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp
index e7fc5c78c0..e24a45b196 100644
--- a/src/nrnoc/treeset.cpp
+++ b/src/nrnoc/treeset.cpp
@@ -1595,7 +1595,6 @@ static double diam_from_list(Section* sec, int inode, Prop* p, double rparent)
 
 void v_setup_vectors(void) {
     int inode, i;
-    NrnThread* _nt;
 
     if (tree_changed) {
         setup_topology(); /* now classical secorder */
@@ -1658,7 +1657,7 @@ void v_setup_vectors(void) {
     reorder_secorder();
 #endif
 
-    FOR_THREADS(_nt) {
+    for (NrnThread* _nt: for_threads(nrn_threads, nrn_nthread)) {
         for (inode = 0; inode < _nt->end; inode++) {
             if (_nt->_v_parent[inode] != NULL) {
                 _nt->_v_parent_index[inode] = _nt->_v_parent[inode]->v_node_index;
@@ -1770,8 +1769,7 @@ void v_setup_vectors(void) {
 
 
 void nrn_matrix_node_free() {
-    NrnThread* nt;
-    FOR_THREADS(nt) {
+    for (NrnThread* nt: for_threads(nrn_threads, nrn_nthread)) {
         if (nt->_sp13_rhs) {
             free(std::exchange(nt->_sp13_rhs, nullptr));
         }
@@ -1785,7 +1783,6 @@ void nrn_matrix_node_free() {
 
 /* 0 means no model, 1 means ODE, 2 means DAE */
 int nrn_modeltype(void) {
-    NrnThread* nt;
     int type;
     v_setup_vectors();
 
@@ -1796,9 +1793,10 @@ int nrn_modeltype(void) {
     type = 0;
     if (nrn_global_ncell > 0) {
         type = 1;
-        FOR_THREADS(nt) if (nt->_ecell_memb_list) {
-            type = 2;
-        }
+        for (const NrnThread* nt: for_threads(nrn_threads, nrn_nthread))
+            if (nt->_ecell_memb_list) {
+                type = 2;
+            }
     }
     if (type == 0 && nrn_nonvint_block_ode_count(0, 0)) {
         type = 1;
@@ -1849,10 +1847,9 @@ and therefore is passed to spSolve as actual_rhs intead of actual_rhs-1.
 
 static void nrn_matrix_node_alloc(void) {
     int i;
-    NrnThread* nt;
 
     nrn_method_consistent();
-    nt = nrn_threads;
+    NrnThread* nt = nrn_threads;
     /*printf("use_sparse13=%d sp13mat=%lx rhs=%lx\n", use_sparse13, (long)nt->_sp13mat,
      * (long)nt->_actual_rhs);*/
     if (use_sparse13) {
@@ -1929,7 +1926,7 @@ static void nrn_matrix_node_alloc(void) {
         }
         nrndae_alloc();
     } else {
-        FOR_THREADS(nt) {
+        for (NrnThread* nt: for_threads(nrn_threads, nrn_nthread)) {
             assert(nrndae_extra_eqn_count() == 0);
             assert(!nt->_ecell_memb_list || nt->_ecell_memb_list->nodecount == 0);
         }
@@ -1972,8 +1969,7 @@ static void nrn_sort_mech_data(
         std::size_t global_i{}, trivial_counter{};
         std::vector<std::size_t> mech_data_permutation(mech_data_size,
                                                        std::numeric_limits<std::size_t>::max());
-        NrnThread* nt{};
-        FOR_THREADS(nt) {
+        for (NrnThread* nt: for_threads(nrn_threads, nrn_nthread)) {
             // the Memb_list for this mechanism in this thread, this might be
             // null if there are no entries, or if it's an artificial cell type(?)
             auto* const ml = nt->_ml_list[type];
@@ -2136,8 +2132,7 @@ static void nrn_sort_node_data(neuron::container::Node::storage::frozen_token_ty
                                                    std::numeric_limits<std::size_t>::max());
     // Process threads one at a time -- this means that the data for each
     // NrnThread will be contiguous.
-    NrnThread* nt{};
-    FOR_THREADS(nt) {
+    for (NrnThread* nt: for_threads(nrn_threads, nrn_nthread)) {
         // What offset in the global node data structure do the values for this thread
         // start at
         nt->_node_data_offset = global_i;
diff --git a/src/nrnpython/convert_cxx_exceptions.hpp b/src/nrnpython/convert_cxx_exceptions.hpp
index 683356e934..a4a8b41f72 100644
--- a/src/nrnpython/convert_cxx_exceptions.hpp
+++ b/src/nrnpython/convert_cxx_exceptions.hpp
@@ -29,7 +29,7 @@ struct error_value_impl<
 
 template <class F, class... Args>
 struct convert_cxx_exceptions_trait {
-    using return_type = typename std::result_of<F(Args...)>::type;
+    using return_type = std::invoke_result_t<F, Args...>;
 
     static return_type error_value() {
         return detail::error_value_impl<return_type>::value();
diff --git a/src/nrnpython/nrnpython.cpp b/src/nrnpython/nrnpython.cpp
index 217679fca6..5dfa174401 100644
--- a/src/nrnpython/nrnpython.cpp
+++ b/src/nrnpython/nrnpython.cpp
@@ -194,9 +194,6 @@ static int nrnpython_start(int b) {
     static int started = 0;
     if (b == 1 && !started) {
         p_nrnpy_pyrun = nrnpy_pyrun;
-        if (nrnpy_nositeflag) {
-            Py_NoSiteFlag = 1;
-        }
         // Create a Python configuration, see
         // https://docs.python.org/3.8/c-api/init_config.html#python-configuration, so that
         // {nrniv,special} -python behaves as similarly as possible to python. In particular this
@@ -204,6 +201,9 @@ static int nrnpython_start(int b) {
         // handle settings like LC_ALL=C, so using a different configuration can lead to surprising
         // differences.
         PythonConfigWrapper config;
+        if (nrnpy_nositeflag) {
+            config->site_import = 0;
+        }
         auto const check = [](const char* desc, PyStatus status) {
             if (PyStatus_Exception(status)) {
                 std::ostringstream oss;