diff --git a/cpp/examples/d_abm.cpp b/cpp/examples/d_abm.cpp
index 09ffa04d5c..4f91f65de0 100644
--- a/cpp/examples/d_abm.cpp
+++ b/cpp/examples/d_abm.cpp
@@ -1,5 +1,5 @@
 /*
-* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC)
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: Julia Bicker, René Schmieding
 *
diff --git a/cpp/examples/graph_abm.cpp b/cpp/examples/graph_abm.cpp
index 3e7443fbc9..dedc28b60e 100644
--- a/cpp/examples/graph_abm.cpp
+++ b/cpp/examples/graph_abm.cpp
@@ -1,5 +1,5 @@
 /*
-* Copyright (C) 2020-2024 MEmilio
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: Julia Bicker
 *
diff --git a/cpp/examples/smm.cpp b/cpp/examples/smm.cpp
index a55785e130..cf13a0c7be 100644
--- a/cpp/examples/smm.cpp
+++ b/cpp/examples/smm.cpp
@@ -1,7 +1,7 @@
 /*
-* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC)
+* Copyright (C) 2020-2025 MEmilio
 *
-* Authors: Julia Bicker, René Schmieding
+* Authors: Julia Bicker, René Schmieding, Kilian Volmer
 *
 * Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
 *
@@ -18,6 +18,7 @@
 * limitations under the License.
 */
 
+#include "memilio/epidemiology/age_group.h"
 #include "smm/simulation.h"
 #include "smm/model.h"
 #include "smm/parameters.h"
@@ -33,44 +34,55 @@ enum class InfectionState
     R,
     D,
     Count
-
 };
 
+using Age     = mio::AgeGroup;
+using Species = mio::AgeGroup;
+
 int main()
 {
 
     //Example how to run the stochastic metapopulation models with four regions
-    const size_t num_regions = 4;
-    using Model              = mio::smm::Model<ScalarType, num_regions, InfectionState>;
+    const size_t num_regions    = 4;
+    const size_t num_age_groups = 1;
+    const size_t num_groups     = 1;
+    using Model                 = mio::smm::Model<ScalarType, num_regions, InfectionState, num_age_groups, num_groups>;
 
     ScalarType numE = 12, numC = 4, numI = 12, numR = 0, numD = 0;
 
     Model model;
     //Population are distributed uniformly to the four regions
     for (size_t r = 0; r < num_regions; ++r) {
-        model.populations[{mio::regions::Region(r), InfectionState::S}] =
+        model.populations[{mio::regions::Region(r), InfectionState::S, Age(0), Species(0)}] =
             (1000 - numE - numC - numI - numR - numD) / num_regions;
-        model.populations[{mio::regions::Region(r), InfectionState::E}] = numE / num_regions;
-        model.populations[{mio::regions::Region(r), InfectionState::C}] = numC / num_regions;
-        model.populations[{mio::regions::Region(r), InfectionState::I}] = numI / num_regions;
-        model.populations[{mio::regions::Region(r), InfectionState::R}] = numR / num_regions;
-        model.populations[{mio::regions::Region(r), InfectionState::D}] = numD / num_regions;
+        model.populations[{mio::regions::Region(r), InfectionState::E, Age(0), Species(0)}] = numE / num_regions;
+        model.populations[{mio::regions::Region(r), InfectionState::C, Age(0), Species(0)}] = numC / num_regions;
+        model.populations[{mio::regions::Region(r), InfectionState::I, Age(0), Species(0)}] = numI / num_regions;
+        model.populations[{mio::regions::Region(r), InfectionState::R, Age(0), Species(0)}] = numR / num_regions;
+        model.populations[{mio::regions::Region(r), InfectionState::D, Age(0), Species(0)}] = numD / num_regions;
     }
 
     //Set infection state adoption and spatial transition rates
-    std::vector<mio::AdoptionRate<ScalarType, InfectionState>> adoption_rates;
-    std::vector<mio::smm::TransitionRate<ScalarType, InfectionState>> transition_rates;
+    std::vector<mio::AdoptionRate<ScalarType, InfectionState, Age, Species>> adoption_rates;
+    std::vector<mio::smm::TransitionRate<ScalarType, InfectionState, Age, Species>> transition_rates;
     for (size_t r = 0; r < num_regions; ++r) {
         adoption_rates.push_back({InfectionState::S,
                                   InfectionState::E,
                                   mio::regions::Region(r),
                                   0.1,
-                                  {{InfectionState::C, 1}, {InfectionState::I, 0.5}}});
-        adoption_rates.push_back({InfectionState::E, InfectionState::C, mio::regions::Region(r), 1.0 / 5., {}});
-        adoption_rates.push_back({InfectionState::C, InfectionState::R, mio::regions::Region(r), 0.2 / 3., {}});
-        adoption_rates.push_back({InfectionState::C, InfectionState::I, mio::regions::Region(r), 0.8 / 3., {}});
-        adoption_rates.push_back({InfectionState::I, InfectionState::R, mio::regions::Region(r), 0.99 / 5., {}});
-        adoption_rates.push_back({InfectionState::I, InfectionState::D, mio::regions::Region(r), 0.01 / 5., {}});
+                                  {{InfectionState::C, 1, mio::regions::Region(3), {Age(0), Species(0)}},
+                                   {InfectionState::I, 0.5, mio::regions::Region(1), {Age(0), Species(0)}}},
+                                  {Age(0), Species(0)}});
+        adoption_rates.push_back(
+            {InfectionState::C, InfectionState::R, mio::regions::Region(r), 0.2 / 3., {}, {Age(0), Species(0)}});
+        adoption_rates.push_back(
+            {InfectionState::E, InfectionState::C, mio::regions::Region(r), 1.0 / 5., {}, {Age(0), Species(0)}});
+        adoption_rates.push_back(
+            {InfectionState::C, InfectionState::I, mio::regions::Region(r), 0.8 / 3., {}, {Age(0), Species(0)}});
+        adoption_rates.push_back(
+            {InfectionState::I, InfectionState::R, mio::regions::Region(r), 0.99 / 5., {}, {Age(0), Species(0)}});
+        adoption_rates.push_back(
+            {InfectionState::I, InfectionState::D, mio::regions::Region(r), 0.01 / 5., {}, {Age(0), Species(0)}});
     }
 
     //Agents in infection state D do not transition
@@ -78,16 +90,24 @@ int main()
         for (size_t i = 0; i < num_regions; ++i) {
             for (size_t j = 0; j < num_regions; ++j)
                 if (i != j) {
-                    transition_rates.push_back(
-                        {InfectionState(s), mio::regions::Region(i), mio::regions::Region(j), 0.01});
-                    transition_rates.push_back(
-                        {InfectionState(s), mio::regions::Region(j), mio::regions::Region(i), 0.01});
+                    transition_rates.push_back({InfectionState(s),
+                                                mio::regions::Region(i),
+                                                mio::regions::Region(j),
+                                                0.01,
+                                                {Age(0), Species(0)},
+                                                {Age(0), Species(0)}});
+                    transition_rates.push_back({InfectionState(s),
+                                                mio::regions::Region(j),
+                                                mio::regions::Region(i),
+                                                0.01,
+                                                {Age(0), Species(0)},
+                                                {Age(0), Species(0)}});
                 }
         }
     }
 
-    model.parameters.get<mio::smm::AdoptionRates<ScalarType, InfectionState>>()   = adoption_rates;
-    model.parameters.get<mio::smm::TransitionRates<ScalarType, InfectionState>>() = transition_rates;
+    model.parameters.get<mio::smm::AdoptionRates<ScalarType, InfectionState, Age, Species>>()   = adoption_rates;
+    model.parameters.get<mio::smm::TransitionRates<ScalarType, InfectionState, Age, Species>>() = transition_rates;
 
     ScalarType dt   = 0.1;
     ScalarType tmax = 30.0;
diff --git a/cpp/memilio/epidemiology/adoption_rate.h b/cpp/memilio/epidemiology/adoption_rate.h
index 705f729025..f1a89a8d2c 100644
--- a/cpp/memilio/epidemiology/adoption_rate.h
+++ b/cpp/memilio/epidemiology/adoption_rate.h
@@ -1,7 +1,7 @@
 /* 
 * Copyright (C) 2020-2025 MEmilio
 *
-* Authors: René Schmieding, Julia Bicker
+* Authors: René Schmieding, Julia Bicker, Kilian Volmer
 *
 * Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
 *
@@ -23,35 +23,45 @@
 #include "memilio/utils/index.h"
 #include "memilio/config.h"
 #include "memilio/geography/regions.h"
+#include <limits>
+#include <tuple>
+#include <optional>
 
 namespace mio
 {
 
-/**
- * @brief Struct defining an influence for a second-order adoption.
- * The population having "status" is multiplied with "factor."
- * @tparam Status An infection state enum.
- */
-template <typename FP, class Status>
-struct Influence {
-    Status status;
-    FP factor;
-};
-
 /**
  * @brief Struct defining a possible status adoption in a Model based on Poisson Processes.
  * The AdoptionRate is considered to be of second-order if there are any "influences".
  * In the d_abm and smm simulations, "from" is implicitly an influence, scaled by "factor". This is multiplied by
  * the sum over all "influences", which scale their "status" with the respective "factor".
  * @tparam Status An infection state enum.
+ * @tparam Groups Additional grouping indices.
  */
-template <typename FP, class Status>
+template <typename FP, class Status, class... Groups>
 struct AdoptionRate {
+
+    /**
+ * @brief Struct defining an influence for a second-order adoption.
+ * The population having "status" is multiplied with "factor."
+ * @tparam status An infection state enum.
+ * @tparam factor Scaling factor for the influence.
+ * @tparam 
+ * @tparam Groups Additional grouping indices.
+ */
+    struct Influence {
+        Status status;
+        FP factor;
+        std::optional<mio::regions::Region> region = std::nullopt;
+        std::tuple<Groups...> group_indices{};
+    };
+
     Status from; // i
     Status to; // j
     mio::regions::Region region; // k
     FP factor; // gammahat_{ij}^k
-    std::vector<Influence<FP, Status>> influences; // influences[tau] = ( Psi_{i,j,tau} , gamma_{i,j,tau} )
+    std::vector<Influence> influences; // influences[tau] = ( Psi_{i,j,tau} , gamma_{i,j,tau} )
+    std::tuple<Groups...> group_indices{};
 };
 
 } // namespace mio
diff --git a/cpp/models/d_abm/model.cpp b/cpp/models/d_abm/model.cpp
index f52bd1c298..528fa1a4f4 100644
--- a/cpp/models/d_abm/model.cpp
+++ b/cpp/models/d_abm/model.cpp
@@ -1,5 +1,5 @@
 /*
-* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC)
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: René Schmieding
 *
diff --git a/cpp/models/d_abm/model.h b/cpp/models/d_abm/model.h
index b0983c5496..5b858a4efa 100644
--- a/cpp/models/d_abm/model.h
+++ b/cpp/models/d_abm/model.h
@@ -1,5 +1,5 @@
 /*
-* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC)
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: René Schmieding, Julia Bicker
 *
diff --git a/cpp/models/d_abm/parameters.h b/cpp/models/d_abm/parameters.h
index bd02befd9b..5b16eb433f 100644
--- a/cpp/models/d_abm/parameters.h
+++ b/cpp/models/d_abm/parameters.h
@@ -1,5 +1,5 @@
 /*
-* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC)
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: René Schmieding, Julia Bicker
 *
diff --git a/cpp/models/d_abm/quad_well.h b/cpp/models/d_abm/quad_well.h
index 305617304e..ef1d0e021e 100644
--- a/cpp/models/d_abm/quad_well.h
+++ b/cpp/models/d_abm/quad_well.h
@@ -1,5 +1,5 @@
 /* 
-* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC)
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: René Schmieding, Julia Bicker
 *
diff --git a/cpp/models/d_abm/simulation.cpp b/cpp/models/d_abm/simulation.cpp
index 3b538490b8..3915b84300 100644
--- a/cpp/models/d_abm/simulation.cpp
+++ b/cpp/models/d_abm/simulation.cpp
@@ -1,5 +1,5 @@
 /*
-* Copyright (C) 2020-2024 MEmilio
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: René Schmieding, Julia Bicker
 *
diff --git a/cpp/models/d_abm/simulation.h b/cpp/models/d_abm/simulation.h
index 54ccad8b65..e20f40f595 100644
--- a/cpp/models/d_abm/simulation.h
+++ b/cpp/models/d_abm/simulation.h
@@ -1,5 +1,5 @@
 /*
-* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC)
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: René Schmieding, Julia Bicker
 *
diff --git a/cpp/models/d_abm/single_well.h b/cpp/models/d_abm/single_well.h
index 62ef7671c6..7bd8cb82e1 100644
--- a/cpp/models/d_abm/single_well.h
+++ b/cpp/models/d_abm/single_well.h
@@ -1,5 +1,5 @@
 /* 
-* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC)
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: Julia Bicker, René Schmieding
 *
diff --git a/cpp/models/graph_abm/graph_abm_mobility.cpp b/cpp/models/graph_abm/graph_abm_mobility.cpp
index af3fae0c89..c8482f2e48 100644
--- a/cpp/models/graph_abm/graph_abm_mobility.cpp
+++ b/cpp/models/graph_abm/graph_abm_mobility.cpp
@@ -1,5 +1,5 @@
 /*
-* Copyright (C) 2020-2024 MEmilio
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: Julia Bicker
 *
diff --git a/cpp/models/graph_abm/graph_abm_mobility.h b/cpp/models/graph_abm/graph_abm_mobility.h
index 2d7a749edf..1af88d311a 100644
--- a/cpp/models/graph_abm/graph_abm_mobility.h
+++ b/cpp/models/graph_abm/graph_abm_mobility.h
@@ -1,5 +1,5 @@
 /*
-* Copyright (C) 2020-2024 MEmilio
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: Julia Bicker
 *
diff --git a/cpp/models/graph_abm/graph_abmodel.h b/cpp/models/graph_abm/graph_abmodel.h
index 0be1a7f222..040cd8c2da 100644
--- a/cpp/models/graph_abm/graph_abmodel.h
+++ b/cpp/models/graph_abm/graph_abmodel.h
@@ -1,5 +1,5 @@
 /*
-* Copyright (C) 2020-2024 MEmilio
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: Julia Bicker
 *
diff --git a/cpp/models/hybrid/temporal_hybrid_model.cpp b/cpp/models/hybrid/temporal_hybrid_model.cpp
index f356bd2321..337b8da3de 100644
--- a/cpp/models/hybrid/temporal_hybrid_model.cpp
+++ b/cpp/models/hybrid/temporal_hybrid_model.cpp
@@ -1,5 +1,5 @@
 /* 
-* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC)
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: Julia Bicker, René Schmieding
 *
diff --git a/cpp/models/hybrid/temporal_hybrid_model.h b/cpp/models/hybrid/temporal_hybrid_model.h
index 25ea35d8d0..3ba096c115 100644
--- a/cpp/models/hybrid/temporal_hybrid_model.h
+++ b/cpp/models/hybrid/temporal_hybrid_model.h
@@ -1,5 +1,5 @@
 /* 
-* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC)
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: Julia Bicker, René Schmieding
 *
diff --git a/cpp/models/smm/model.cpp b/cpp/models/smm/model.cpp
index af0d58eb73..d613003ac2 100644
--- a/cpp/models/smm/model.cpp
+++ b/cpp/models/smm/model.cpp
@@ -1,5 +1,5 @@
 /*
-* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC)
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: René Schmieding
 *
diff --git a/cpp/models/smm/model.h b/cpp/models/smm/model.h
index a30a6f365e..b5942bc251 100644
--- a/cpp/models/smm/model.h
+++ b/cpp/models/smm/model.h
@@ -1,7 +1,7 @@
 /*
-* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC)
+* Copyright (C) 2020-2025 MEmilio
 *
-* Authors: René Schmieding, Julia Bicker
+* Authors: René Schmieding, Julia Bicker, Kilian Volmer
 *
 * Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
 *
@@ -22,6 +22,7 @@
 #define MIO_SMM_MODEL_H
 
 #include "memilio/config.h"
+#include "memilio/epidemiology/age_group.h"
 #include "smm/parameters.h"
 #include "memilio/compartments/compartmental_model.h"
 #include "memilio/epidemiology/populations.h"
@@ -32,21 +33,29 @@ namespace mio
 namespace smm
 {
 
+template <class T, std::size_t>
+using age_group = T;
+
 /**
  * @brief Stochastic Metapopulation Model.
  * @tparam regions Number of regions.
  * @tparam Status An infection state enum.
  */
-template <typename FP, size_t regions, class Status>
-class Model : public mio::CompartmentalModel<FP, Status, mio::Populations<FP, mio::regions::Region, Status>,
-                                             ParametersBase<FP, Status>>
+template <typename FP, size_t regions, class Status, size_t... Groups>
+class Model : public mio::CompartmentalModel<
+                  FP, Status, mio::Populations<FP, mio::regions::Region, Status, age_group<mio::AgeGroup, Groups>...>,
+                  ParametersBase<FP, Status, age_group<mio::AgeGroup, Groups>...>>
 {
-    using Base = mio::CompartmentalModel<FP, Status, mio::Populations<FP, mio::regions::Region, Status>,
-                                         ParametersBase<FP, Status>>;
+    using Base =
+        mio::CompartmentalModel<FP, Status,
+                                mio::Populations<FP, mio::regions::Region, Status, age_group<mio::AgeGroup, Groups>...>,
+                                ParametersBase<FP, Status, age_group<mio::AgeGroup, Groups>...>>;
+    using Index = mio::Index<mio::regions::Region, Status, age_group<mio::AgeGroup, Groups>...>;
 
 public:
     Model()
-        : Base(typename Base::Populations({static_cast<mio::regions::Region>(regions), Status::Count}, 0.0),
+        : Base(typename Base::Populations(
+                   {static_cast<mio::regions::Region>(regions), Status::Count, static_cast<mio::AgeGroup>(Groups)...}),
                typename Base::ParameterSet())
     {
     }
@@ -57,26 +66,46 @@ class Model : public mio::CompartmentalModel<FP, Status, mio::Populations<FP, mi
      * @param[in] x The current state of the model.
      * @return Current value of the adoption rate.
      */
-    FP evaluate(const AdoptionRate<FP, Status>& rate, const Eigen::VectorX<FP>& x) const
+    FP evaluate(const AdoptionRate<FP, Status, age_group<mio::AgeGroup, Groups>...>& rate,
+                const Eigen::VectorXd& x) const
     {
-        const auto& pop   = this->populations;
-        const auto source = pop.get_flat_index({rate.region, rate.from});
+        const auto& pop       = this->populations;
+        const auto index_from = std::apply(
+            [&](auto&&... args) {
+                return Index{rate.region, rate.from, std::forward<decltype(args)>(args)...};
+            },
+            rate.group_indices);
+        const auto source = pop.get_flat_index(index_from); // Why is here rate.from used? KV
         // determine order and calculate rate
         if (rate.influences.size() == 0) { // first order adoption
             return rate.factor * x[source];
         }
         else { // second order adoption
-            FP N = 0;
-            for (size_t s = 0; s < static_cast<size_t>(Status::Count); ++s) {
-                N += x[pop.get_flat_index({rate.region, Status(s)})];
-            }
             // accumulate influences
             FP influences = 0.0;
             for (size_t i = 0; i < rate.influences.size(); i++) {
-                influences +=
-                    rate.influences[i].factor * x[pop.get_flat_index({rate.region, rate.influences[i].status})];
+                FP N = 0; // Welches N brauchen wir hier??
+
+                for (size_t s = 0; s < static_cast<size_t>(Status::Count); ++s) {
+                    const auto index = std::apply(
+                        [&](auto&&... args) {
+                            return Index{rate.influences[i].region.value_or(rate.region), Status(s),
+                                         std::forward<decltype(args)>(args)...};
+                        },
+                        rate.influences[i].group_indices);
+                    N += x[pop.get_flat_index(index)];
+                }
+                const auto index = std::apply(
+                    [&](auto&&... args) {
+                        return Index{rate.influences[i].region.value_or(rate.region), rate.influences[i].status,
+                                     std::forward<decltype(args)>(args)...};
+                    },
+                    rate.influences[i].group_indices);
+                if (N > 0) {
+                    influences += rate.influences[i].factor * x[pop.get_flat_index(index)] / N;
+                }
             }
-            return (N > 0) ? (rate.factor * x[source] * influences / N) : 0;
+            return rate.factor * x[source] * influences;
         }
     }
 
@@ -86,9 +115,15 @@ class Model : public mio::CompartmentalModel<FP, Status, mio::Populations<FP, mi
      * @param[in] x The current state of the model.
      * @return Current value of the transition rate.
      */
-    FP evaluate(const TransitionRate<FP, Status>& rate, const Eigen::VectorX<FP>& x) const
+    FP evaluate(const TransitionRate<FP, Status, age_group<mio::AgeGroup, Groups>...>& rate,
+                const Eigen::VectorXd& x) const
     {
-        const auto source = this->populations.get_flat_index({rate.from, rate.status});
+        auto index = std::apply(
+            [&](auto&&... args) {
+                return Index{rate.from, rate.status, std::forward<decltype(args)>(args)...};
+            },
+            rate.group_indices_from);
+        const auto source = this->populations.get_flat_index(index);
         return rate.factor * x[source];
     }
 
diff --git a/cpp/models/smm/parameters.h b/cpp/models/smm/parameters.h
index 4cc54ebb5b..fe27c35526 100644
--- a/cpp/models/smm/parameters.h
+++ b/cpp/models/smm/parameters.h
@@ -1,5 +1,5 @@
 /*
-* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC)
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: René Schmieding, Julia Bicker
 *
@@ -25,6 +25,7 @@
 #include "memilio/geography/regions.h"
 #include "memilio/utils/parameter_set.h"
 #include "memilio/epidemiology/adoption_rate.h"
+#include <tuple>
 
 namespace mio
 {
@@ -35,9 +36,9 @@ namespace smm
  * @brief A vector of AdoptionRate%s, see mio::AdoptionRate
  * @tparam Status An infection state enum.
  */
-template <typename FP, class Status>
+template <typename FP, class Status, class... Groups>
 struct AdoptionRates {
-    using Type = std::vector<AdoptionRate<FP, Status>>;
+    using Type = std::vector<AdoptionRate<FP, Status, Groups...>>;
     const static std::string name()
     {
         return "AdoptionRates";
@@ -48,25 +49,26 @@ struct AdoptionRates {
  * @brief Struct defining a possible regional transition in a Model based on Poisson Processes.
  * @tparam Status An infection state enum.
  */
-template <typename FP, class Status>
+template <typename FP, class Status, class... Groups>
 struct TransitionRate {
     Status status; // i
     mio::regions::Region from; // k
     mio::regions::Region to; // l
     FP factor; // lambda_i^{kl}
+    std::tuple<Groups...> group_indices_from{};
+    std::tuple<Groups...> group_indices_to{};
 };
-
-template <typename FP, class Status>
+template <typename FP, class Status, class... Groups>
 struct TransitionRates {
-    using Type = std::vector<TransitionRate<FP, Status>>;
+    using Type = std::vector<TransitionRate<FP, Status, Groups...>>;
     const static std::string name()
     {
         return "TransitionRates";
     }
 };
 
-template <typename FP, class Status>
-using ParametersBase = mio::ParameterSet<AdoptionRates<FP, Status>, TransitionRates<FP, Status>>;
+template <typename FP, class Status, class... Groups>
+using ParametersBase = mio::ParameterSet<AdoptionRates<FP, Status, Groups...>, TransitionRates<FP, Status, Groups...>>;
 
 } // namespace smm
 
diff --git a/cpp/models/smm/simulation.h b/cpp/models/smm/simulation.h
index 4231f76a05..8fb1423571 100644
--- a/cpp/models/smm/simulation.h
+++ b/cpp/models/smm/simulation.h
@@ -1,5 +1,5 @@
 /*
-* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC)
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: René Schmieding, Julia Bicker
 *
@@ -22,9 +22,13 @@
 #define MIO_SMM_SIMULATION_H
 
 #include "memilio/config.h"
+#include "memilio/epidemiology/age_group.h"
+#include "memilio/utils/index.h"
+#include "memilio/utils/logging.h"
 #include "smm/model.h"
 #include "smm/parameters.h"
 #include "memilio/compartments/simulation.h"
+#include <utility>
 
 namespace mio
 {
@@ -37,12 +41,12 @@ namespace smm
  * @tparam regions The number of regions.
  * @tparam Status An infection state enum.
  */
-template <typename FP, size_t regions, class Status>
+template <typename FP, size_t regions, class Status, size_t... Groups>
 class Simulation
 {
 public:
-public:
-    using Model = smm::Model<FP, regions, Status>;
+    using Model = smm::Model<FP, regions, Status, Groups...>;
+    using Index = mio::Index<mio::regions::Region, Status, age_group<mio::AgeGroup, Groups>...>;
 
     /**
      * @brief Set up the simulation for a Stochastic Metapopulation Model.
@@ -100,18 +104,38 @@ class Simulation
             if (next_event < adoption_rates().size()) {
                 // perform adoption event
                 const auto& rate = adoption_rates()[next_event];
-                m_result.get_last_value()[m_model->populations.get_flat_index({rate.region, rate.from})] -= 1;
-                m_model->populations[{rate.region, rate.from}] -= 1.0;
-                m_result.get_last_value()[m_model->populations.get_flat_index({rate.region, rate.to})] += 1;
-                m_model->populations[{rate.region, rate.to}] += 1.0;
+                auto index_from  = std::apply(
+                    [&](auto&&... args) {
+                        return Index{rate.region, rate.from, std::forward<decltype(args)>(args)...};
+                    },
+                    rate.group_indices);
+                m_result.get_last_value()[m_model->populations.get_flat_index(index_from)] -= 1;
+                m_model->populations[index_from] -= 1;
+                auto index_to = std::apply(
+                    [&](auto&&... args) {
+                        return Index{rate.region, rate.to, std::forward<decltype(args)>(args)...};
+                    },
+                    rate.group_indices);
+                m_result.get_last_value()[m_model->populations.get_flat_index(index_to)] += 1;
+                m_model->populations[index_to] += 1;
             }
             else {
                 // perform transition event
                 const auto& rate = transition_rates()[next_event - adoption_rates().size()];
-                m_result.get_last_value()[m_model->populations.get_flat_index({rate.from, rate.status})] -= 1;
-                m_model->populations[{rate.from, rate.status}] -= 1.0;
-                m_result.get_last_value()[m_model->populations.get_flat_index({rate.to, rate.status})] += 1;
-                m_model->populations[{rate.to, rate.status}] += 1.0;
+                auto index_from  = std::apply(
+                    [&](auto&&... args) {
+                        return Index{rate.from, rate.status, std::forward<decltype(args)>(args)...};
+                    },
+                    rate.group_indices_from);
+                m_result.get_last_value()[m_model->populations.get_flat_index(index_from)] -= 1;
+                m_model->populations[index_from] -= 1;
+                auto index_to = std::apply(
+                    [&](auto&&... args) {
+                        return Index{rate.to, rate.status, std::forward<decltype(args)>(args)...};
+                    },
+                    rate.group_indices_to);
+                m_result.get_last_value()[m_model->populations.get_flat_index(index_to)] += 1;
+                m_model->populations[index_to] += 1;
             }
             // update internal times
             for (size_t i = 0; i < m_internal_time.size(); i++) {
@@ -129,6 +153,7 @@ class Simulation
             m_result.get_last_value() = m_result[m_result.get_num_time_points() - 2];
         }
         return m_result.get_last_value();
+        // }
     }
 
     /**
@@ -160,21 +185,25 @@ class Simulation
     /**
      * @brief Returns the model's transition rates.
      */
-    inline constexpr const typename smm::TransitionRates<FP, Status>::Type& transition_rates()
+    inline constexpr const typename smm::TransitionRates<FP, Status, age_group<mio::AgeGroup, Groups>...>::Type&
+    transition_rates()
     {
-        return m_model->parameters.template get<smm::TransitionRates<FP, Status>>();
+        return m_model->parameters
+            .template get<smm::TransitionRates<FP, Status, age_group<mio::AgeGroup, Groups>...>>();
     }
 
     /**
      * @brief Returns the model's adoption rates.
      */
-    inline constexpr const typename smm::AdoptionRates<FP, Status>::Type& adoption_rates()
+    inline constexpr const typename smm::AdoptionRates<FP, Status, age_group<mio::AgeGroup, Groups>...>::Type&
+    adoption_rates()
     {
-        return m_model->parameters.template get<smm::AdoptionRates<FP, Status>>();
+        return m_model->parameters.template get<smm::AdoptionRates<FP, Status, age_group<mio::AgeGroup, Groups>...>>();
     }
 
     /**
-     * @brief Calculate current values for m_current_rates and m_waiting_times.
+     * @brief 
+     * 
      */
     inline void update_current_rates_and_waiting_times()
     {
@@ -213,7 +242,7 @@ class Simulation
     std::vector<FP> m_current_rates; ///< Current values of both types of rates i.e. adoption and transition rates.
 };
 
-} //namespace smm
+} // namespace smm
 } // namespace mio
 
 #endif
diff --git a/cpp/tests/test_d_abm_model.cpp b/cpp/tests/test_d_abm_model.cpp
index c31b4aa5e4..79397d8e65 100644
--- a/cpp/tests/test_d_abm_model.cpp
+++ b/cpp/tests/test_d_abm_model.cpp
@@ -1,5 +1,5 @@
 /* 
-* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC)
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: Julia Bicker
 *
diff --git a/cpp/tests/test_smm_model.cpp b/cpp/tests/test_smm_model.cpp
index 30ba70e0b5..9e7be226f4 100644
--- a/cpp/tests/test_smm_model.cpp
+++ b/cpp/tests/test_smm_model.cpp
@@ -1,7 +1,7 @@
 /*
-* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC)
+* Copyright (C) 2020-2025 MEmilio
 *
-* Authors: Julia Bicker
+* Authors: Julia Bicker, Kilian Volmer
 *
 * Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
 *
@@ -23,7 +23,9 @@
 #include <cstddef>
 #include <cstdint>
 #include <gtest/gtest.h>
+#include "memilio/epidemiology/age_group.h"
 #include "memilio/utils/compiler_diagnostics.h"
+#include "memilio/utils/logging.h"
 #include "smm/model.h"
 #include "smm/parameters.h"
 #include "smm/simulation.h"
@@ -45,6 +47,14 @@ enum class InfectionState
 
 };
 
+enum class Infections
+{
+    S,
+    I,
+    R,
+    Count
+};
+
 TEST(TestSMM, evaluateAdoptionRate)
 {
     //Test whether the adoption rates are evaluated correctly.
@@ -80,31 +90,80 @@ TEST(TestSMM, evaluateAdoptionRate)
     EXPECT_EQ(model.evaluate(adoption_rates[1], model.populations.get_compartments()), 2.);
 }
 
+TEST(TestSMM, evaluateinterregionalAdoptionRate)
+{
+    using Model = mio::smm::Model<double, 2, Infections, 2>;
+    using Age   = mio::AgeGroup;
+
+    std::vector<mio::AdoptionRate<double, Infections, Age>> adoption_rates;
+    //Second-order adoption
+    adoption_rates.push_back({Infections::S,
+                              Infections::I,
+                              mio::regions::Region(0),
+                              0.1,
+                              {{Infections::I, 0.1, mio::regions::Region(0), {Age(1)}},
+                               {Infections::I, 0.2, mio::regions::Region(1), {Age(1)}}},
+                              Age(0)});
+    //First-order adoption
+    adoption_rates.push_back({Infections::I, Infections::R, mio::regions::Region(0), 0.2, {}, {Age(1)}});
+    Model model;
+
+    model.populations[{mio::regions::Region(0), Infections::S, Age(0)}] = 50;
+    model.populations[{mio::regions::Region(0), Infections::I, Age(0)}] = 10;
+    model.populations[{mio::regions::Region(0), Infections::R, Age(0)}] = 0;
+
+    model.populations[{mio::regions::Region(0), Infections::S, Age(1)}] = 100;
+    model.populations[{mio::regions::Region(0), Infections::I, Age(1)}] = 20;
+    model.populations[{mio::regions::Region(0), Infections::R, Age(1)}] = 0;
+
+    model.populations[{mio::regions::Region(1), Infections::S, Age(0)}] = 40;
+    model.populations[{mio::regions::Region(1), Infections::I, Age(0)}] = 80;
+    model.populations[{mio::regions::Region(1), Infections::R, Age(0)}] = 0;
+
+    model.populations[{mio::regions::Region(1), Infections::S, Age(1)}] = 80;
+    model.populations[{mio::regions::Region(1), Infections::I, Age(1)}] = 16;
+    model.populations[{mio::regions::Region(1), Infections::R, Age(1)}] = 0;
+
+    EXPECT_DOUBLE_EQ(model.evaluate(adoption_rates[0], model.populations.get_compartments()),
+                     0.1 * 50. * (0.1 * 20. * 1. / 120. + 0.2 * 16 * 1. / 96.));
+    EXPECT_DOUBLE_EQ(model.evaluate(adoption_rates[1], model.populations.get_compartments()), 4.);
+}
+
 TEST(TestSMM, evaluateTransitionRate)
 {
     //Same test as 'evaluateAdoptionRate' only for spatial transition rates.
     //Transition rates are given by: rate.factor * N(rate.status, rate.from)
-    using Model = mio::smm::Model<double, 2, InfectionState>;
+    using Model = mio::smm::Model<double, 2, InfectionState, 2>;
 
     Model model;
     //Initialize model populations
-    model.populations[{mio::regions::Region(0), InfectionState::S}] = 50;
-    model.populations[{mio::regions::Region(0), InfectionState::E}] = 10;
-    model.populations[{mio::regions::Region(0), InfectionState::C}] = 5;
-    model.populations[{mio::regions::Region(0), InfectionState::I}] = 0;
-    model.populations[{mio::regions::Region(0), InfectionState::R}] = 0;
-    model.populations[{mio::regions::Region(0), InfectionState::D}] = 0;
-
-    model.populations[{mio::regions::Region(1), InfectionState::S}] = 55;
-    model.populations[{mio::regions::Region(1), InfectionState::E}] = 10;
-    model.populations[{mio::regions::Region(1), InfectionState::C}] = 0;
-    model.populations[{mio::regions::Region(1), InfectionState::I}] = 0;
-    model.populations[{mio::regions::Region(1), InfectionState::R}] = 0;
-    model.populations[{mio::regions::Region(1), InfectionState::D}] = 0;
+    model.populations[{mio::regions::Region(0), InfectionState::S, mio::AgeGroup(0)}] = 50;
+    model.populations[{mio::regions::Region(0), InfectionState::E, mio::AgeGroup(0)}] = 10;
+    model.populations[{mio::regions::Region(0), InfectionState::C, mio::AgeGroup(0)}] = 5;
+    model.populations[{mio::regions::Region(0), InfectionState::I, mio::AgeGroup(0)}] = 0;
+    model.populations[{mio::regions::Region(0), InfectionState::R, mio::AgeGroup(0)}] = 0;
+    model.populations[{mio::regions::Region(0), InfectionState::D, mio::AgeGroup(0)}] = 0;
+
+    model.populations[{mio::regions::Region(1), InfectionState::S, mio::AgeGroup(0)}] = 55;
+    model.populations[{mio::regions::Region(1), InfectionState::E, mio::AgeGroup(0)}] = 10;
+    model.populations[{mio::regions::Region(1), InfectionState::C, mio::AgeGroup(0)}] = 0;
+    model.populations[{mio::regions::Region(1), InfectionState::I, mio::AgeGroup(0)}] = 0;
+    model.populations[{mio::regions::Region(1), InfectionState::R, mio::AgeGroup(0)}] = 0;
+    model.populations[{mio::regions::Region(1), InfectionState::D, mio::AgeGroup(0)}] = 0;
     //Set transition rates
-    std::vector<mio::smm::TransitionRate<double, InfectionState>> transition_rates;
-    transition_rates.push_back({InfectionState::S, mio::regions::Region(0), mio::regions::Region(1), 0.01});
-    transition_rates.push_back({InfectionState::E, mio::regions::Region(1), mio::regions::Region(0), 0.1});
+    std::vector<mio::smm::TransitionRate<double, InfectionState, mio::AgeGroup>> transition_rates;
+    transition_rates.push_back({InfectionState::S,
+                                mio::regions::Region(0),
+                                mio::regions::Region(1),
+                                0.01,
+                                {mio::AgeGroup(0)},
+                                {mio::AgeGroup(0)}});
+    transition_rates.push_back({InfectionState::E,
+                                mio::regions::Region(1),
+                                mio::regions::Region(0),
+                                0.1,
+                                {mio::AgeGroup(0)},
+                                {mio::AgeGroup(0)}});
 
     EXPECT_EQ(model.evaluate(transition_rates[0], model.populations.get_compartments()), 0.5);
     EXPECT_EQ(model.evaluate(transition_rates[1], model.populations.get_compartments()), 1.);
diff --git a/pycode/memilio-epidata/memilio/epidata/getJHData.py b/pycode/memilio-epidata/memilio/epidata/getJHData.py
index 19ec2b67f9..dd0a3bb873 100644
--- a/pycode/memilio-epidata/memilio/epidata/getJHData.py
+++ b/pycode/memilio-epidata/memilio/epidata/getJHData.py
@@ -1,5 +1,5 @@
 #############################################################################
-# Copyright (C) 2020-2024 MEmilio
+# Copyright (C) 2020-2025 MEmilio
 #
 # Authors: Kathrin Rack
 #
diff --git a/pycode/memilio-generation/memilio/generation/template/template_cpp.txt b/pycode/memilio-generation/memilio/generation/template/template_cpp.txt
index 28d9f80042..e33668d1ee 100644
--- a/pycode/memilio-generation/memilio/generation/template/template_cpp.txt
+++ b/pycode/memilio-generation/memilio/generation/template/template_cpp.txt
@@ -1,5 +1,5 @@
 /* 
-* Copyright (C) 2020-2024 MEmilio
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: Martin Siggel, Daniel Abele, Martin J. Kuehn, Jan Kleinert
 *
diff --git a/pycode/memilio-generation/memilio/generation/template/template_py.txt b/pycode/memilio-generation/memilio/generation/template/template_py.txt
index 2ae2f78e5a..d1cc9710c7 100644
--- a/pycode/memilio-generation/memilio/generation/template/template_py.txt
+++ b/pycode/memilio-generation/memilio/generation/template/template_py.txt
@@ -1,5 +1,5 @@
 #############################################################################
-# Copyright (C) 2020-2024 MEmilio
+# Copyright (C) 2020-2025 MEmilio
 #
 # Authors: Daniel Abele
 #
diff --git a/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.cpp.txt b/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.cpp.txt
index b419020ba3..edacafb623 100644
--- a/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.cpp.txt
+++ b/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.cpp.txt
@@ -1,5 +1,5 @@
 /* 
-* Copyright (C) 2020-2024 MEmilio
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: Martin Siggel, Daniel Abele, Martin J. Kuehn, Jan Kleinert
 *
diff --git a/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.py.txt b/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.py.txt
index 3f289a1251..1fd1861e46 100644
--- a/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.py.txt
+++ b/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.py.txt
@@ -1,5 +1,5 @@
 #############################################################################
-# Copyright (C) 2020-2024 MEmilio
+# Copyright (C) 2020-2025 MEmilio
 #
 # Authors: Daniel Abele
 #