SciCompMod · HenrZu · Oct 27, 2025 · Aug 14, 2025 · Sep 5, 2025 · Sep 8, 2025
diff --git a/cpp/examples/sde_sirs.cpp b/cpp/examples/sde_sirs.cpp
@@ -48,7 +48,8 @@ int main()
     model.parameters.get<mio::ssirs::ContactPatterns<ScalarType>>().get_baseline()(0, 0) = 20.7;
     model.parameters.get<mio::ssirs::ContactPatterns<ScalarType>>().add_damping(0.6,
                                                                                 mio::SimulationTime<ScalarType>(12.5));
-
+    model.parameters.set<mio::ssirs::StartDay<ScalarType>>(60);
+    model.parameters.set<mio::ssirs::Seasonality<ScalarType>>(0.2);
     model.check_constraints();
 
     auto ssirs = mio::simulate_stochastic(t0, tmax, dt, model);

diff --git a/cpp/models/sde_sirs/model.h b/cpp/models/sde_sirs/model.h
@@ -56,7 +56,12 @@ class Model
                    Eigen::Ref<Eigen::VectorX<FP>> flows) const
     {
         auto& params = Base::parameters;
-        FP coeffStoI = params.template get<ContactPatterns<FP>>().get_matrix_at(SimulationTime<FP>(t))(0, 0) *
+        // effective contact rate by contact rate between groups i and j and damping j
+        FP season_val =
+            (1 + params.template get<Seasonality<FP>>() *
+                     sin(std::numbers::pi_v<ScalarType> * ((params.template get<StartDay<FP>>() + t) / 182.5 + 0.5)));
+
+        FP coeffStoI = season_val * params.template get<ContactPatterns<FP>>().get_matrix_at(SimulationTime<FP>(t))(0, 0) *
                        params.template get<TransmissionProbabilityOnContact<FP>>() / Base::populations.get_total();
 
         flows[this->template get_flat_flow_index<InfectionState::Susceptible, InfectionState::Infected>()] =

diff --git a/cpp/models/sde_sirs/parameters.h b/cpp/models/sde_sirs/parameters.h
@@ -98,9 +98,46 @@ struct ContactPatterns {
     }
 };
 
+/**
+ * @brief The start day in the SIRS model
+ * The start day defines in which season the simulation can be started
+ * If the start day is 180 and simulation takes place from t0=0 to
+ * tmax=100 the days 180 to 280 of the year are simulated
+ */
+template <typename FP>
+struct StartDay {
+    using Type = FP;
+    static Type get_default()
+    {
+        return Type(0.0);
+    }
+    static std::string name()
+    {
+        return "StartDay";
+    }
+};
+
+/**
+ * @brief The seasonality in the SIRS model.
+ * The seasonality is given as (1+k*sin()) where the sine
+ * curve is below one in summer and above one in winter
+ */
+template <typename FP>
+struct Seasonality {
+    using Type = UncertainValue<FP>;
+    static Type get_default()
+    {
+        return Type(0.);
+    }
+    static std::string name()
+    {
+        return "Seasonality";
+    }
+};
+
 template <typename FP>
 using ParametersBase =
-    ParameterSet<TransmissionProbabilityOnContact<FP>, TimeInfected<FP>, ContactPatterns<FP>, TimeImmune<FP>>;
+    ParameterSet<TransmissionProbabilityOnContact<FP>, TimeInfected<FP>, ContactPatterns<FP>, TimeImmune<FP>, Seasonality<FP>, StartDay<FP>>;
 
 /**
  * @brief Parameters of SIR model.
@@ -132,6 +169,12 @@ class Parameters : public ParametersBase<FP>
         FP tol_times = 1e-1;
 
         int corrected = false;
+        if (this->template get<Seasonality<FP>>() < 0.0 || this->template get<Seasonality<FP>>() > 0.5) {
+            log_warning("Constraint check: Parameter Seasonality changed from {:0.4f} to {:d}",
+                        this->template get<Seasonality<FP>>(), 0);
+            this->template set<Seasonality<FP>>(0);
+            corrected = true;
+        }
         if (this->template get<TimeInfected<FP>>() < tol_times) {
             log_warning("Constraint check: Parameter TimeInfected changed from {:.4f} to {:.4f}. Please note that "
                         "unreasonably small compartment stays lead to massively increased run time. Consider to cancel "
@@ -167,6 +210,10 @@ class Parameters : public ParametersBase<FP>
     {
         FP tol_times = 1e-1;
 
+        if (this->template get<Seasonality<FP>>() < 0.0 || this->template get<Seasonality<FP>>() > 0.5) {
+            log_error("Constraint check: Parameter Seasonality smaller {:d} or larger {:d}", 0, 0.5);
+            return true;
+        }
         if (this->template get<TimeInfected<FP>>() < tol_times) {
             log_error("Constraint check: Parameter TimeInfected {:.4f} smaller or equal {:.4f}. Please note that "
                       "unreasonably small compartment stays lead to massively increased run time. Consider to cancel "

diff --git a/cpp/tests/test_sde_sirs.cpp b/cpp/tests/test_sde_sirs.cpp
@@ -90,6 +90,8 @@ TEST(TestSdeSirs, check_constraints_parameters)
     parameters.set<mio::ssirs::TimeImmune<double>>(6);
     parameters.set<mio::ssirs::TransmissionProbabilityOnContact<double>>(0.04);
     parameters.get<mio::ssirs::ContactPatterns<double>>().get_baseline()(0, 0) = 10;
+    parameters.set<mio::ssirs::StartDay<double>>(30);
+    parameters.set<mio::ssirs::Seasonality<double>>(0.3);
 
     // model.check_constraints() combines the functions from population and parameters.
     // We only want to test the functions for the parameters defined in parameters.h
@@ -107,6 +109,10 @@ TEST(TestSdeSirs, check_constraints_parameters)
     parameters.set<mio::ssirs::TimeImmune<double>>(6);
     parameters.set<mio::ssirs::TransmissionProbabilityOnContact<double>>(10.);
     EXPECT_EQ(parameters.check_constraints(), 1);
+
+    parameters.set<mio::ssirs::TransmissionProbabilityOnContact<double>>(0.04);
+    parameters.set<mio::ssirs::Seasonality<double>>(-2.);
+    EXPECT_EQ(parameters.check_constraints(), 1);
     mio::set_log_level(mio::LogLevel::warn);
 }
 
@@ -119,6 +125,8 @@ TEST(TestSdeSirs, apply_constraints_parameters)
     parameters.set<mio::ssirs::TimeImmune<double>>(6);
     parameters.set<mio::ssirs::TransmissionProbabilityOnContact<double>>(0.04);
     parameters.get<mio::ssirs::ContactPatterns<double>>().get_baseline()(0, 0) = 10;
+    parameters.set<mio::ssirs::StartDay<double>>(30);
+    parameters.set<mio::ssirs::Seasonality<double>>(0.3);
 
     EXPECT_EQ(parameters.apply_constraints(), 0);
 
@@ -135,5 +143,9 @@ TEST(TestSdeSirs, apply_constraints_parameters)
     parameters.set<mio::ssirs::TransmissionProbabilityOnContact<double>>(10.);
     EXPECT_EQ(parameters.apply_constraints(), 1);
     EXPECT_NEAR(parameters.get<mio::ssirs::TransmissionProbabilityOnContact<double>>(), 0.0, 1e-14);
+
+    parameters.set<mio::ssirs::Seasonality<double>>(-2.);
+    EXPECT_EQ(parameters.apply_constraints(), 1);
+    EXPECT_NEAR(parameters.get<mio::ssirs::Seasonality<double>>(), 0.0, 1e-14);
     mio::set_log_level(mio::LogLevel::warn);
 }
diff --git a/docs/source/cpp/models/ssirs.rst b/docs/source/cpp/models/ssirs.rst
@@ -16,19 +16,21 @@ Below is an overview of the model architecture and its compartments.
 .. image:: https://martinkuehn.eu/research/images/sirs.png
    :alt: SIR_model
 
-+-------------------------------+-----------------------------------------------+--------------------------------------------------------------------------------------------------+
-| Mathematical variable         | C++ variable name                             | Description                                                                                      |
-+===============================+===============================================+==================================================================================================+
-| :math:`\phi`                  | ``ContactPatterns``                           | Daily contact rate / Number of daily contacts.                                                   |
-+-------------------------------+-----------------------------------------------+--------------------------------------------------------------------------------------------------+
-| :math:`\rho`                  | ``TransmissionProbabilityOnContact``          | Transmission risk for people located in the Susceptible compartment.                             |
-+-------------------------------+-----------------------------------------------+--------------------------------------------------------------------------------------------------+
-| :math:`N`                     | ``populations.get_total()``                   | Total population.                                                                                |
-+-------------------------------+-----------------------------------------------+--------------------------------------------------------------------------------------------------+
-| :math:`T_{I}`                 | ``TimeInfected``                              | Time in days an individual stays in the Infected compartment.                                    |
-+-------------------------------+-----------------------------------------------+--------------------------------------------------------------------------------------------------+
-| :math:`T_{R}`                 | ``TimeImmune``                                | Time in days an individual stays in the Recovered compartment before becoming Susceptible again. |
-+-------------------------------+-----------------------------------------------+--------------------------------------------------------------------------------------------------+
++-------------------------------+-----------------------------------------------+------------------------------------------------------------------------------------------------------------+
+| Mathematical variable         | C++ variable name                             | Description                                                                                                |
++===============================+===============================================+============================================================================================================+
+| :math:`\phi`                  | ``ContactPatterns``                           | Daily contact rate / Number of daily contacts.                                                             |
++-------------------------------+-----------------------------------------------+------------------------------------------------------------------------------------------------------------+
+| :math:`\rho`                  | ``TransmissionProbabilityOnContact``          | Transmission risk for people located in the Susceptible compartment.                                       |
++-------------------------------+-----------------------------------------------+------------------------------------------------------------------------------------------------------------+
+| :math:`N`                     | ``populations.get_total()``                   | Total population.                                                                                          |
++-------------------------------+-----------------------------------------------+------------------------------------------------------------------------------------------------------------+
+| :math:`T_{I}`                 | ``TimeInfected``                              | Time in days an individual stays in the Infected compartment.                                              |
++-------------------------------+-----------------------------------------------+------------------------------------------------------------------------------------------------------------+
+| :math:`T_{R}`                 | ``TimeImmune``                                | Time in days an individual stays in the Recovered compartment before becoming Susceptible again.           |
++-------------------------------+-----------------------------------------------+------------------------------------------------------------------------------------------------------------+
+| :math:`k`                     | ``Seasonality``                               | Influence of the seasons is given by :math:`s_k(t) = 1 + k \sin \left(\frac{t}{182.5} + \frac{1}{2}\right)`|
++-------------------------------+-----------------------------------------------+------------------------------------------------------------------------------------------------------------+
 
 An example can be found in the
 `examples/ode_sir.cpp <https://github.com/SciCompMod/memilio/blob/main/cpp/examples/sde_sirs.cpp>`_.

diff --git a/pycode/examples/simulation/sde_sir_simple.py b/pycode/examples/simulation/sde_sir_simple.py
@@ -0,0 +1,71 @@
+#############################################################################
+# Copyright (C) 2020-2025 MEmilio
+#
+# Authors: Maximilian Betz
+#
+# Contact: Martin J. Kuehn <[email protected]>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#############################################################################
+import argparse
+
+import numpy as np
+
+from memilio.simulation import AgeGroup, Damping
+from memilio.simulation.ssir import InfectionState as State
+from memilio.simulation.ssir import (Model, simulate_stochastic)
+
+
+def run_sde_sir_simulation():
+    """Runs SDE SIR model"""
+
+    tmax = 5.  # simulation time frame
+    dt = 0.1
+
+    # Initialize model
+    model = Model()
+
+    # Mean time in Infected compartment
+    model.parameters.TimeInfected.value = 10.
+
+    model.parameters.TransmissionProbabilityOnContact.value = 1.
+
+    # Initial number of people per compartment
+    total_population = 10000
+    model.populations[State.Infected] = 100
+    model.populations[State.Recovered] = 1000
+    model.populations.set_difference_from_total(
+        (State.Susceptible), total_population)
+
+    model.parameters.ContactPatterns.baseline = np.ones(
+        (1, 1)) * 2.7
+    model.parameters.ContactPatterns.minimum = np.zeros(
+        (1, 1))
+    model.parameters.ContactPatterns.add_damping(
+        Damping(coeffs=np.r_[0.6], t=2., level=0, type=0))
+
+    # Check parameter constraints
+    model.check_constraints()
+
+    # Run Simulation
+    result = simulate_stochastic(0., tmax, dt, model)
+
+    result.print_table(False, ["S", "I", "R"], 16, 5)
+
+
+if __name__ == "__main__":
+    arg_parser = argparse.ArgumentParser(
+        'sde_sir_simple',
+        description='Simple example demonstrating the setup and simulation of the SDE SIR model.')
+    args = arg_parser.parse_args()
+    run_sde_sir_simulation()
diff --git a/pycode/examples/simulation/sde_sirs_simple.py b/pycode/examples/simulation/sde_sirs_simple.py
@@ -0,0 +1,76 @@
+#############################################################################
+# Copyright (C) 2020-2025 MEmilio
+#
+# Authors: Maximilian Betz
+#
+# Contact: Martin J. Kuehn <[email protected]>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#############################################################################
+import argparse
+
+import numpy as np
+
+from memilio.simulation import AgeGroup, Damping
+from memilio.simulation.ssirs import InfectionState as State
+from memilio.simulation.ssirs import (
+    Model, simulate_stochastic, interpolate_simulation_result)
+
+
+def run_sde_sirs_simulation():
+    """Runs SDE SIRS model"""
+
+    tmax = 5.  # simulation time frame
+    dt = 0.001
+
+    # Initialize Model
+    model = Model()
+
+    # Mean time in Infected compartment
+    model.parameters.TimeInfected.value = 10.
+    model.parameters.TimeImmune.value = 100.
+
+    model.parameters.TransmissionProbabilityOnContact.value = 1.
+
+    # Initial number of people per compartment
+    total_population = 10000
+    model.populations[State.Infected] = 100
+    model.populations[State.Recovered] = 1000
+    model.populations.set_difference_from_total(
+        (State.Susceptible), total_population)
+
+    model.parameters.ContactPatterns.baseline = np.ones(
+        (1, 1)) * 20.7
+    model.parameters.ContactPatterns.minimum = np.zeros(
+        (1, 1))
+    model.parameters.ContactPatterns.add_damping(
+        Damping(coeffs=np.r_[0.6], t=2, level=0, type=0))
+
+    # Check parameter constraints
+    model.check_constraints()
+
+    # Run Simulation
+    result = simulate_stochastic(0., days, dt, model)
+
+    # Interpolate results
+    result = interpolate_simulation_result(result)
+
+    result.print_table(False, ["Susceptible", "Infected", "Recovered"], 16, 5)
+
+
+if __name__ == "__main__":
+    arg_parser = argparse.ArgumentParser(
+        'sde_sirs_simple',
+        description='Simple example demonstrating the setup and simulation of the SDE SIRS model.')
+    args = arg_parser.parse_args()
+    run_sde_sirs_simulation()
diff --git a/pycode/memilio-simulation/CMakeLists.txt b/pycode/memilio-simulation/CMakeLists.txt
@@ -87,6 +87,15 @@ add_pymio_module(_simulation_osecirvvs
     SOURCES memilio/simulation/bindings/models/osecirvvs.cpp
 )
 
+add_pymio_module(_simulation_ssir
+    LINKED_LIBRARIES memilio sde_sir
+    SOURCES memilio/simulation/bindings/models/ssir.cpp
+)
+
+add_pymio_module(_simulation_ssirs
+    LINKED_LIBRARIES memilio sde_sirs
+    SOURCES memilio/simulation/bindings/models/ssirs.cpp
+)
 add_pymio_module(_simulation_omseirs4
     LINKED_LIBRARIES memilio ode_mseirs4
     SOURCES memilio/simulation/bindings/models/omseirs4.cpp

diff --git a/pycode/memilio-simulation/memilio/simulation/__init__.py b/pycode/memilio-simulation/memilio/simulation/__init__.py
@@ -45,6 +45,13 @@ def __getattr__(attr):
     elif attr == "osecirvvs":
         import memilio.simulation.osecirvvs as osecirvvs
         return osecirvvs
+    elif attr == "ssir":
+        import memilio.simulation.ssir as ssir
+        return ssir
+
+    elif attr == "ssirs":
+        import memilio.simulation.ssirs as ssirs
+        return ssirs
 
     raise AttributeError("module {!r} has no attribute "
                          "{!r}".format(__name__, attr))