diff --git a/examples/ActivityDrivenInertial/conf.toml b/examples/ActivityDrivenInertial/conf.toml
new file mode 100644
index 0000000..383ba76
--- /dev/null
+++ b/examples/ActivityDrivenInertial/conf.toml
@@ -0,0 +1,34 @@
+[simulation]
+model = "ActivityDrivenInertial"
+# rng_seed = 120 # Leaving this empty will pick a random seed
+
+[io]
+n_output_network = 20 # Write the network every 20 iterations
+n_output_agents = 1   # Write the opinions of agents after every iteration
+print_progress = true # Print the iteration time ; if not set, then does not print
+
+[model]
+max_iterations = 500 # If not set, max iterations is infinite
+
+[ActivityDrivenInertial]
+dt = 0.01               # Timestep for the integration of the coupled ODEs
+m = 10                  # Number of agents contacted, when the agent is active
+eps = 0.01              # Minimum activity epsilon; a_i belongs to [epsilon,1]
+gamma = 2.1             # Exponent of activity power law distribution of activities
+reciprocity = 0.65      # probability that when agent i contacts j via weighted reservoir sampling, j also sends feedback to i. So every agent can have more than m incoming connections
+homophily = 1.0         # aka beta. if zero, agents pick their interaction partners at random
+alpha = 3.0             # Controversialness of the issue, must be greater than 0.
+K = 2.0
+mean_activities = false # Use the mean value of the powerlaw distribution for the activities of all agents
+mean_weights = false    # Use the meanfield approximation of the network edges
+
+reluctances = true       # Assigns a "reluctance" (m_i) to each agent. By default; false and every agent has a reluctance of 1
+reluctance_mean = 1.0    # Mean of distribution before drawing from a truncated normal distribution (default set to 1.0)
+reluctance_sigma = 0.25  # Width of normal distribution (before truncating)
+reluctance_eps = 0.01    # Minimum such that the normal distribution is truncated at this value
+
+friction_coefficient = 1.0 # Friction coefficient, making agents tend to go to rest without acceleration
+
+[network]
+number_of_agents = 1000
+connections_per_agent = 10
diff --git a/include/agents/inertial_agent.hpp b/include/agents/inertial_agent.hpp
new file mode 100644
index 0000000..469376a
--- /dev/null
+++ b/include/agents/inertial_agent.hpp
@@ -0,0 +1,68 @@
+#pragma once
+
+#include "agent.hpp"
+#include "agent_io.hpp"
+#include <util/misc.hpp>
+
+namespace Seldon
+{
+
+struct InertialAgentData
+{
+    double opinion    = 0;   // x_i
+    double activity   = 0;   // alpha_i
+    double reluctance = 1.0; // m_i
+    double velocity   = 0.0; // d(x_i)/dt
+};
+
+using InertialAgent = Agent<InertialAgentData>;
+
+template<>
+inline std::string agent_to_string<InertialAgent>( const InertialAgent & agent )
+{
+    return fmt::format(
+        "{}, {}, {}, {}", agent.data.opinion, agent.data.velocity, agent.data.activity, agent.data.reluctance );
+}
+
+template<>
+inline std::string opinion_to_string<InertialAgent>( const InertialAgent & agent )
+{
+    return fmt::format( "{}", agent.data.opinion );
+}
+
+template<>
+inline InertialAgent agent_from_string<InertialAgent>( const std::string & str )
+{
+    InertialAgent res{};
+
+    auto callback = [&]( int idx_list, std::string & substr )
+    {
+        if( idx_list == 0 )
+        {
+            res.data.opinion = std::stod( substr );
+        }
+        else if( idx_list == 1 )
+        {
+            res.data.velocity = std::stod( substr );
+        }
+        else if( idx_list == 2 )
+        {
+            res.data.activity = std::stod( substr );
+        }
+        else if( idx_list == 3 )
+        {
+            res.data.reluctance = std::stod( substr );
+        }
+    };
+
+    Seldon::parse_comma_separated_list( str, callback );
+
+    return res;
+};
+
+template<>
+inline std::vector<std::string> agent_to_string_column_names<InertialAgent>()
+{
+    return { "opinion", "velocity", "activity", "reluctance" };
+}
+} // namespace Seldon
\ No newline at end of file
diff --git a/include/config_parser.hpp b/include/config_parser.hpp
index 208e27d..326dc59 100644
--- a/include/config_parser.hpp
+++ b/include/config_parser.hpp
@@ -24,7 +24,8 @@ namespace Seldon::Config
 enum class Model
 {
     DeGroot,
-    ActivityDrivenModel,
+    ActivityDrivenModel, // @TODO : no need for model here
+    ActivityDrivenInertial,
     DeffuantModel
 };
 
@@ -81,6 +82,11 @@ struct ActivityDrivenSettings
     double covariance_factor          = 0.0;
 };
 
+struct ActivityDrivenInertialSettings : public ActivityDrivenSettings
+{
+    double friction_coefficient = 1.0;
+};
+
 struct InitialNetworkSettings
 {
     std::optional<std::string> file;
@@ -90,7 +96,8 @@ struct InitialNetworkSettings
 
 struct SimulationOptions
 {
-    using ModelVariantT = std::variant<DeGrootSettings, ActivityDrivenSettings, DeffuantSettings>;
+    using ModelVariantT
+        = std::variant<DeGrootSettings, ActivityDrivenSettings, ActivityDrivenInertialSettings, DeffuantSettings>;
     Model model;
     std::string model_string;
     int rng_seed = std::random_device()();
diff --git a/include/connectivity.hpp b/include/connectivity.hpp
index a9d0159..1af9c8d 100644
--- a/include/connectivity.hpp
+++ b/include/connectivity.hpp
@@ -75,7 +75,7 @@ class TarjanConnectivityAlgo
                 {
                     lowest[v] = std::min( lowest[v], num[u] );
                 } // u not processed
-            } // u has been visited
+            }     // u has been visited
         }
 
         // Now v has been processed
diff --git a/include/model_factory.hpp b/include/model_factory.hpp
index d5864ac..a87f9dd 100644
--- a/include/model_factory.hpp
+++ b/include/model_factory.hpp
@@ -3,6 +3,7 @@
 #include "models/ActivityDrivenModel.hpp"
 #include "models/DeGroot.hpp"
 #include "models/DeffuantModel.hpp"
+#include "models/InertialModel.hpp"
 #include "network.hpp"
 #include <memory>
 #include <random>
@@ -61,6 +62,23 @@ create_model_activity_driven( Network<AgentT> & network, const ModelVariantT & m
     }
 }
 
+template<typename AgentT>
+inline auto create_model_activity_driven_inertial(
+    Network<AgentT> & network, const ModelVariantT & model_settings, std::mt19937 & gen )
+{
+    if constexpr( std::is_same_v<AgentT, InertialModel::AgentT> )
+    {
+        auto settings = std::get<Config::ActivityDrivenInertialSettings>( model_settings );
+        auto model    = std::make_unique<InertialModel>( settings, network, gen );
+        return model;
+    }
+    else
+    {
+        throw std::runtime_error( "Incompatible agent and model type!" );
+        return std::unique_ptr<Model<AgentT>>{};
+    }
+}
+
 template<typename AgentT>
 inline auto create_model_deffuant( Network<AgentT> & network, const ModelVariantT & model_settings, std::mt19937 & gen )
 {
diff --git a/include/models/ActivityDrivenModel.hpp b/include/models/ActivityDrivenModel.hpp
index 06d830f..1b0ab83 100644
--- a/include/models/ActivityDrivenModel.hpp
+++ b/include/models/ActivityDrivenModel.hpp
@@ -1,10 +1,11 @@
 #pragma once
 
 #include "agents/activity_agent.hpp"
+#include "agents/inertial_agent.hpp"
 #include "config_parser.hpp"
-
 #include "model.hpp"
 #include "network.hpp"
+#include "network_generation.hpp"
 #include <cstddef>
 #include <random>
 #include <set>
@@ -15,23 +16,62 @@
 namespace Seldon
 {
 
-class ActivityDrivenModel : public Model<ActivityAgent>
+template<typename AgentT_>
+class ActivityDrivenModelAbstract : public Model<AgentT_>
 {
 public:
-    using AgentT   = ActivityAgent;
+    using AgentT   = AgentT_;
     using NetworkT = Network<AgentT>;
+    using WeightT  = typename NetworkT::WeightT;
 
-    ActivityDrivenModel( const Config::ActivityDrivenSettings & settings, NetworkT & network, std::mt19937 & gen );
+    ActivityDrivenModelAbstract(
+        const Config::ActivityDrivenSettings & settings, NetworkT & network, std::mt19937 & gen )
+            : Model<AgentT>( settings.max_iterations ),
+              network( network ),
+              contact_prob_list( std::vector<std::vector<WeightT>>( network.n_agents() ) ),
+              gen( gen ),
+              dt( settings.dt ),
+              m( settings.m ),
+              eps( settings.eps ),
+              gamma( settings.gamma ),
+              alpha( settings.alpha ),
+              homophily( settings.homophily ),
+              reciprocity( settings.reciprocity ),
+              K( settings.K ),
+              mean_activities( settings.mean_activities ),
+              mean_weights( settings.mean_weights ),
+              use_reluctances( settings.use_reluctances ),
+              reluctance_mean( settings.reluctance_mean ),
+              reluctance_sigma( settings.reluctance_sigma ),
+              reluctance_eps( settings.reluctance_eps ),
+              n_bots( settings.n_bots ),
+              bot_m( settings.bot_m ),
+              bot_activity( settings.bot_activity ),
+              bot_opinion( settings.bot_opinion ),
+              bot_homophily( settings.bot_homophily )
+    {
+        get_agents_from_power_law();
 
-    void iteration() override;
+        if( mean_weights )
+        {
+            auto agents_copy = network.agents;
+            network          = NetworkGeneration::generate_fully_connected<AgentT>( network.n_agents() );
+            network.agents   = agents_copy;
+        }
+    }
 
-private:
+    void iteration() override{};
+
+protected:
     NetworkT & network;
-    std::vector<std::vector<NetworkT::WeightT>> contact_prob_list; // Probability of choosing i in 1 to m rounds
+
+private:
+    std::vector<std::vector<WeightT>> contact_prob_list; // Probability of choosing i in 1 to m rounds
     // Random number generation
     std::mt19937 & gen; // reference to simulation Mersenne-Twister engine
     std::set<std::pair<size_t, size_t>> reciprocal_edge_buffer{};
 
+protected:
     // Model-specific parameters
     double dt{}; // Timestep for the integration of the coupled ODEs
     // Various free parameters
@@ -66,13 +106,222 @@ class ActivityDrivenModel : public Model<ActivityAgent>
     std::vector<double> k3_buffer{};
     std::vector<double> k4_buffer{};
 
-    void get_agents_from_power_law();
+private:
+    void get_agents_from_power_law()
+    {
+        std::uniform_real_distribution<> dis_opinion( -1, 1 ); // Opinion initial values
+        power_law_distribution<> dist_activity( eps, gamma );
+        truncated_normal_distribution<> dist_reluctance( reluctance_mean, reluctance_sigma, reluctance_eps );
+
+        auto mean_activity = dist_activity.mean();
 
+        // Initial conditions for the opinions, initialize to [-1,1]
+        // The activities should be drawn from a power law distribution
+        for( size_t i = 0; i < network.agents.size(); i++ )
+        {
+            network.agents[i].data.opinion = dis_opinion( gen ); // Draw the opinion value
+
+            if( !mean_activities )
+            { // Draw from a power law distribution (1-gamma)/(1-eps^(1-gamma)) * a^(-gamma)
+                network.agents[i].data.activity = dist_activity( gen );
+            }
+            else
+            {
+                network.agents[i].data.activity = mean_activity;
+            }
+
+            if( use_reluctances )
+            {
+                network.agents[i].data.reluctance = dist_reluctance( gen );
+            }
+        }
+
+        if( bot_present() )
+        {
+            for( size_t bot_idx = 0; bot_idx < n_bots; bot_idx++ )
+            {
+                network.agents[bot_idx].data.opinion  = bot_opinion[bot_idx];
+                network.agents[bot_idx].data.activity = bot_activity[bot_idx];
+            }
+        }
+    }
+
+    // The weight for contact between two agents
+    double homophily_weight( size_t idx_contacter, size_t idx_contacted )
+    {
+        double homophily = this->homophily;
+
+        if( idx_contacted == idx_contacter )
+            return 0.0;
+
+        if( bot_present() && idx_contacter < n_bots )
+            homophily = this->bot_homophily[idx_contacter];
+
+        constexpr double tolerance = 1e-10;
+        auto opinion_diff
+            = std::abs( network.agents[idx_contacter].data.opinion - network.agents[idx_contacted].data.opinion );
+        opinion_diff = std::max( tolerance, opinion_diff );
+
+        return std::pow( opinion_diff, -homophily );
+    }
+
+    void update_network_probabilistic()
+    {
+        network.switch_direction_flag();
+
+        std::uniform_real_distribution<> dis_activation( 0.0, 1.0 );
+        std::uniform_real_distribution<> dis_reciprocation( 0.0, 1.0 );
+        std::vector<size_t> contacted_agents{};
+        reciprocal_edge_buffer.clear(); // Clear the reciprocal edge buffer
+        for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
+        {
+            // Test if the agent is activated
+            bool activated = dis_activation( gen ) < network.agents[idx_agent].data.activity;
+
+            if( activated )
+            {
+                // Implement the weight for the probability of agent `idx_agent` contacting agent `j`
+                // Not normalised since this is taken care of by the reservoir sampling
+
+                int m_temp = this->m;
+
+                if( bot_present() && idx_agent < n_bots )
+                {
+                    m_temp = bot_m[idx_agent];
+                }
+
+                reservoir_sampling_A_ExpJ(
+                    m_temp, network.n_agents(), [&]( int j ) { return homophily_weight( idx_agent, j ); },
+                    contacted_agents, gen );
+
+                // Fill the outgoing edges into the reciprocal edge buffer
+                for( const auto & idx_outgoing : contacted_agents )
+                {
+                    reciprocal_edge_buffer.insert(
+                        { idx_agent, idx_outgoing } ); // insert the edge idx_agent -> idx_outgoing
+                }
+
+                // Set the *outgoing* edges
+                network.set_neighbours_and_weights( idx_agent, contacted_agents, 1.0 );
+            }
+            else
+            {
+                network.set_neighbours_and_weights( idx_agent, {}, {} );
+            }
+        }
+
+        // Reciprocity check
+        for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
+        {
+            // Get the outgoing edges
+            auto contacted_agents = network.get_neighbours( idx_agent );
+            // For each outgoing edge we check if the reverse edge already exists
+            for( const auto & idx_outgoing : contacted_agents )
+            {
+                // If the edge is not reciprocated
+                if( !reciprocal_edge_buffer.contains( { idx_outgoing, idx_agent } ) )
+                {
+                    if( dis_reciprocation( gen ) < reciprocity )
+                    {
+                        network.push_back_neighbour_and_weight( idx_outgoing, idx_agent, 1.0 );
+                    }
+                }
+            }
+        }
+
+        network.toggle_incoming_outgoing(); // switch direction, so that we have incoming edges
+    }
+
+    void update_network_mean()
+    {
+        std::vector<WeightT> weights( network.n_agents(), 0.0 );
+
+        // Set all weights to zero in the beginning
+        for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
+        {
+            network.set_weights( idx_agent, weights );
+            contact_prob_list[idx_agent] = weights; // set to zero
+        }
+
+        auto probability_helper = []( double omega, size_t m )
+        {
+            double p = 0;
+            for( size_t i = 1; i <= m; i++ )
+                p += ( std::pow( -omega, i + 1 ) + omega ) / ( omega + 1 );
+            return p;
+        };
+
+        for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
+        {
+            // Implement the weight for the probability of agent `idx_agent` contacting agent `j`
+            // Not normalised since this is taken care of by the reservoir sampling
+
+            double normalization = 0;
+            for( size_t k = 0; k < network.n_agents(); k++ )
+            {
+                normalization += homophily_weight( idx_agent, k );
+            }
+
+            // Go through all the neighbours of idx_agent
+            // Calculate the probability of i contacting j (in 1 to m rounds, assuming
+            // the agent is activated
+            int m_temp = m;
+            if( bot_present() && idx_agent < n_bots )
+            {
+                m_temp = bot_m[idx_agent];
+            }
+
+            double activity = std::max( 1.0, network.agents[idx_agent].data.activity );
+            for( size_t j = 0; j < network.n_agents(); j++ )
+            {
+                double omega = homophily_weight( idx_agent, j ) / normalization;
+                // We can calculate the probability of i contacting j ( i->j )
+                // Update contact prob_list (outgoing)
+                contact_prob_list[idx_agent][j] = activity * probability_helper( omega, m_temp );
+            }
+        }
+
+        for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
+        {
+            // Calculate the actual weights and reciprocity
+            for( size_t j = 0; j < network.n_agents(); j++ )
+            {
+                double prob_contact_ij = contact_prob_list[idx_agent][j]; // outgoing probabilites
+                double prob_contact_ji = contact_prob_list[j][idx_agent];
+
+                // Set the incoming agent weight, j-i in weight list
+                double & win_ji = network.get_weights( j )[idx_agent];
+                win_ji += prob_contact_ij;
+
+                // Handle the reciprocity for j->i
+                // Update incoming weight i-j
+                double & win_ij = network.get_weights( idx_agent )[j];
+
+                // The probability of reciprocating is
+                win_ij += ( 1.0 - prob_contact_ji ) * reciprocity * prob_contact_ij;
+            }
+        }
+    }
+
+protected:
     [[nodiscard]] bool bot_present() const
     {
         return n_bots > 0;
     }
 
+    void update_network()
+    {
+
+        if( !mean_weights )
+        {
+            update_network_probabilistic();
+        }
+        else
+        {
+            update_network_mean();
+        }
+    }
+
     template<typename Opinion_Callback>
     void get_euler_slopes( std::vector<double> & k_buffer, Opinion_Callback opinion )
     {
@@ -93,16 +342,12 @@ class ActivityDrivenModel : public Model<ActivityAgent>
                 k_buffer[idx_agent] += 1.0 / network.agents[idx_agent].data.reluctance * K * weight_buffer[j]
                                        * std::tanh( alpha * opinion( j_index ) );
             }
-            // Multiply by the timestep
-            k_buffer[idx_agent] *= dt;
+            // Here, we won't multiply by the timestep.
+            // Instead multiply in the update rule
         }
     }
-
-    // The weight for contact between two agents
-    double homophily_weight( size_t idx_contacter, size_t idx_contacted );
-    void update_network_probabilistic();
-    void update_network_mean();
-    void update_network();
 };
 
+using ActivityDrivenModel = ActivityDrivenModelAbstract<ActivityAgent>;
+
 } // namespace Seldon
\ No newline at end of file
diff --git a/include/models/InertialModel.hpp b/include/models/InertialModel.hpp
new file mode 100644
index 0000000..ba95a0b
--- /dev/null
+++ b/include/models/InertialModel.hpp
@@ -0,0 +1,40 @@
+#pragma once
+
+#include "agents/activity_agent.hpp"
+#include "agents/inertial_agent.hpp"
+#include "config_parser.hpp"
+#include "model.hpp"
+#include "models/ActivityDrivenModel.hpp"
+#include "network.hpp"
+#include "network_generation.hpp"
+#include <cstddef>
+#include <random>
+#include <set>
+#include <string>
+#include <utility>
+#include <vector>
+
+namespace Seldon
+{
+
+class InertialModel : public ActivityDrivenModelAbstract<InertialAgent>
+{
+public:
+    using AgentT   = InertialAgent;
+    using NetworkT = Network<AgentT>;
+    using WeightT  = typename NetworkT::WeightT;
+
+    InertialModel( const Config::ActivityDrivenInertialSettings & settings, NetworkT & network, std::mt19937 & gen );
+
+    void iteration() override;
+
+private:
+    double friction_coefficient = 1.0;
+    std::vector<double> drift_t_buffer{};
+    std::vector<double> drift_next_t_buffer{};
+
+    void calc_velocity();
+    void calc_position();
+};
+
+} // namespace Seldon
\ No newline at end of file
diff --git a/include/simulation.hpp b/include/simulation.hpp
index 4eca2fc..118b824 100644
--- a/include/simulation.hpp
+++ b/include/simulation.hpp
@@ -70,6 +70,10 @@ class Simulation : public SimulationInterface
         {
             model = ModelFactory::create_model_activity_driven( network, options.model_settings, gen );
         }
+        else if( options.model == Config::Model::ActivityDrivenInertial )
+        {
+            model = ModelFactory::create_model_activity_driven_inertial( network, options.model_settings, gen );
+        }
         else if( options.model == Config::Model::DeffuantModel )
         {
             auto deffuant_settings = std::get<Config::DeffuantSettings>( options.model_settings );
diff --git a/meson.build b/meson.build
index e7017c3..d5a4548 100644
--- a/meson.build
+++ b/meson.build
@@ -12,6 +12,7 @@ sources_seldon = [
   'src/models/ActivityDrivenModel.cpp',
   'src/models/DeffuantModel.cpp',
   'src/models/DeffuantModelVector.cpp',
+  'src/models/InertialModel.cpp',
   'src/util/tomlplusplus.cpp',
 ]
 
@@ -26,6 +27,7 @@ tests = [
   ['Test_Tarjan', 'test/test_tarjan.cpp'],
   ['Test_DeGroot', 'test/test_deGroot.cpp'],
   ['Test_Activity_Driven', 'test/test_activity.cpp'],
+  ['Test_Activity_Driven_Inertial', 'test/test_activity_inertial.cpp'],
   ['Test_Deffuant', 'test/test_deffuant.cpp'],
   ['Test_Network', 'test/test_network.cpp'],
   ['Test_Network_Generation', 'test/test_network_generation.cpp'],
diff --git a/src/config_parser.cpp b/src/config_parser.cpp
index b65d83a..e0e054e 100644
--- a/src/config_parser.cpp
+++ b/src/config_parser.cpp
@@ -29,9 +29,82 @@ Model model_string_to_enum( std::string_view model_string )
     {
         return Model::DeffuantModel;
     }
+    else if( model_string == "ActivityDrivenInertial" )
+    {
+        return Model::ActivityDrivenInertial;
+    }
     throw std::runtime_error( fmt::format( "Invalid model string {}", model_string ) );
 }
 
+void set_if_specified( auto & opt, const auto & toml_opt )
+{
+    using T    = typename std::remove_reference<decltype( opt )>::type;
+    auto t_opt = toml_opt.template value<T>();
+    if( t_opt.has_value() )
+        opt = t_opt.value();
+}
+
+// Convenience function to parse activity settings
+void parse_activity_settings( auto & model_settings, const auto & toml_model_opt, const auto & tbl )
+{
+    set_if_specified( model_settings.dt, toml_model_opt["dt"] );
+    set_if_specified( model_settings.m, toml_model_opt["m"] );
+    set_if_specified( model_settings.eps, toml_model_opt["eps"] );
+    set_if_specified( model_settings.gamma, toml_model_opt["gamma"] );
+    set_if_specified( model_settings.homophily, toml_model_opt["homophily"] );
+    set_if_specified( model_settings.reciprocity, toml_model_opt["reciprocity"] );
+    set_if_specified( model_settings.alpha, toml_model_opt["alpha"] );
+    set_if_specified( model_settings.K, toml_model_opt["K"] );
+    // Mean activity model options
+    set_if_specified( model_settings.mean_activities, toml_model_opt["mean_activities"] );
+    set_if_specified( model_settings.mean_weights, toml_model_opt["mean_weights"] );
+    // Reluctances
+    set_if_specified( model_settings.use_reluctances, toml_model_opt["reluctances"] );
+    set_if_specified( model_settings.reluctance_mean, toml_model_opt["reluctance_mean"] );
+    set_if_specified( model_settings.reluctance_sigma, toml_model_opt["reluctance_sigma"] );
+    set_if_specified( model_settings.reluctance_eps, toml_model_opt["reluctance_eps"] );
+
+    model_settings.max_iterations = tbl["model"]["max_iterations"].template value<int>();
+
+    // bot
+    set_if_specified( model_settings.n_bots, toml_model_opt["n_bots"] );
+
+    auto push_back_bot_array = [&]( auto toml_node, auto & options_array, auto default_value )
+    {
+        if( toml_node.is_array() )
+        {
+            toml::array * toml_arr = toml_node.as_array();
+
+            toml_arr->for_each(
+                [&]( auto && elem )
+                {
+                    if( elem.is_integer() )
+                    {
+                        options_array.push_back( elem.as_integer()->get() );
+                    }
+                    else if( elem.is_floating_point() )
+                    {
+                        options_array.push_back( elem.as_floating_point()->get() );
+                    }
+                } );
+        }
+        else
+        {
+            options_array = std::vector<decltype( default_value )>( model_settings.n_bots, default_value );
+        }
+    };
+
+    auto bot_opinion   = toml_model_opt["bot_opinion"];
+    auto bot_m         = toml_model_opt["bot_m"];
+    auto bot_activity  = toml_model_opt["bot_activity"];
+    auto bot_homophily = toml_model_opt["bot_homophily"];
+
+    push_back_bot_array( bot_m, model_settings.bot_m, model_settings.m );
+    push_back_bot_array( bot_opinion, model_settings.bot_opinion, 0.0 );
+    push_back_bot_array( bot_activity, model_settings.bot_activity, 1.0 );
+    push_back_bot_array( bot_homophily, model_settings.bot_homophily, model_settings.homophily );
+}
+
 SimulationOptions parse_config_file( std::string_view config_file_path )
 {
     SimulationOptions options;
@@ -41,14 +114,6 @@ SimulationOptions parse_config_file( std::string_view config_file_path )
 
     options.rng_seed = tbl["simulation"]["rng_seed"].value_or( int( options.rng_seed ) );
 
-    auto set_if_specified = [&]( auto & opt, const auto & toml_opt )
-    {
-        using T    = typename std::remove_reference<decltype( opt )>::type;
-        auto t_opt = toml_opt.template value<T>();
-        if( t_opt.has_value() )
-            opt = t_opt.value();
-    };
-
     // Parse output settings
     options.output_settings.n_output_network = tbl["io"]["n_output_network"].value<size_t>();
     options.output_settings.n_output_agents  = tbl["io"]["n_output_agents"].value<size_t>();
@@ -87,63 +152,15 @@ SimulationOptions parse_config_file( std::string_view config_file_path )
     {
         auto model_settings = ActivityDrivenSettings();
 
-        set_if_specified( model_settings.dt, tbl[options.model_string]["dt"] );
-        set_if_specified( model_settings.m, tbl[options.model_string]["m"] );
-        set_if_specified( model_settings.eps, tbl[options.model_string]["eps"] );
-        set_if_specified( model_settings.gamma, tbl[options.model_string]["gamma"] );
-        set_if_specified( model_settings.homophily, tbl[options.model_string]["homophily"] );
-        set_if_specified( model_settings.reciprocity, tbl[options.model_string]["reciprocity"] );
-        set_if_specified( model_settings.alpha, tbl[options.model_string]["alpha"] );
-        set_if_specified( model_settings.K, tbl[options.model_string]["K"] );
-        // Mean activity model options
-        set_if_specified( model_settings.mean_activities, tbl[options.model_string]["mean_activities"] );
-        set_if_specified( model_settings.mean_weights, tbl[options.model_string]["mean_weights"] );
-        // Reluctances
-        set_if_specified( model_settings.use_reluctances, tbl[options.model_string]["reluctances"] );
-        set_if_specified( model_settings.reluctance_mean, tbl[options.model_string]["reluctance_mean"] );
-        set_if_specified( model_settings.reluctance_sigma, tbl[options.model_string]["reluctance_sigma"] );
-        set_if_specified( model_settings.reluctance_eps, tbl[options.model_string]["reluctance_eps"] );
-
-        model_settings.max_iterations = tbl["model"]["max_iterations"].value<int>();
-
-        // bot
-        set_if_specified( model_settings.n_bots, tbl[options.model_string]["n_bots"] );
-
-        auto push_back_bot_array = [&]( auto toml_node, auto & options_array, auto default_value )
-        {
-            if( toml_node.is_array() )
-            {
-                toml::array * toml_arr = toml_node.as_array();
-
-                toml_arr->for_each(
-                    [&]( auto && elem )
-                    {
-                        if( elem.is_integer() )
-                        {
-                            options_array.push_back( elem.as_integer()->get() );
-                        }
-                        else if( elem.is_floating_point() )
-                        {
-                            options_array.push_back( elem.as_floating_point()->get() );
-                        }
-                    } );
-            }
-            else
-            {
-                options_array = std::vector<decltype( default_value )>( model_settings.n_bots, default_value );
-            }
-        };
-
-        auto bot_opinion   = tbl[options.model_string]["bot_opinion"];
-        auto bot_m         = tbl[options.model_string]["bot_m"];
-        auto bot_activity  = tbl[options.model_string]["bot_activity"];
-        auto bot_homophily = tbl[options.model_string]["bot_homophily"];
-
-        push_back_bot_array( bot_m, model_settings.bot_m, model_settings.m );
-        push_back_bot_array( bot_opinion, model_settings.bot_opinion, 0.0 );
-        push_back_bot_array( bot_activity, model_settings.bot_activity, 1.0 );
-        push_back_bot_array( bot_homophily, model_settings.bot_homophily, model_settings.homophily );
+        parse_activity_settings( model_settings, tbl[options.model_string], tbl );
+        options.model_settings = model_settings;
+    }
+    else if( options.model == Model::ActivityDrivenInertial )
+    {
+        auto model_settings = ActivityDrivenInertialSettings();
 
+        parse_activity_settings( model_settings, tbl[options.model_string], tbl );
+        set_if_specified( model_settings.friction_coefficient, tbl[options.model_string]["friction_coefficient"] );
         options.model_settings = model_settings;
     }
 
@@ -183,10 +200,8 @@ void validate_settings( const SimulationOptions & options )
     // @TODO: Check that start_output is less than the max_iterations?
     check( name_and_var( options.output_settings.start_output ), g_zero );
 
-    if( options.model == Model::ActivityDrivenModel )
+    auto validate_activity = [&]( const auto & model_settings )
     {
-        auto model_settings = std::get<ActivityDrivenSettings>( options.model_settings );
-
         check( name_and_var( model_settings.dt ), g_zero );
         check( name_and_var( model_settings.m ), geq_zero );
         check( name_and_var( model_settings.eps ), g_zero );
@@ -207,6 +222,19 @@ void validate_settings( const SimulationOptions & options )
         check( name_and_var( model_settings.bot_activity ), check_bot_size, bot_msg );
         check( name_and_var( model_settings.bot_opinion ), check_bot_size, bot_msg );
         check( name_and_var( model_settings.bot_homophily ), check_bot_size, bot_msg );
+    };
+
+    if( options.model == Model::ActivityDrivenModel )
+    {
+        auto model_settings = std::get<ActivityDrivenSettings>( options.model_settings );
+
+        validate_activity( model_settings );
+    }
+    else if( options.model == Model::ActivityDrivenInertial )
+    {
+        auto model_settings = std::get<ActivityDrivenInertialSettings>( options.model_settings );
+        check( name_and_var( model_settings.friction_coefficient ), geq_zero );
+        validate_activity( model_settings );
     }
     else if( options.model == Model::DeGroot )
     {
@@ -225,7 +253,8 @@ void validate_settings( const SimulationOptions & options )
         {
             const std::string basic_deff_msg
                 = "The basic Deffuant model has not been implemented with multiple dimensions";
-            check( name_and_var( model_settings.dim ), []( auto x ) { return x == 1; }, basic_deff_msg );
+            check(
+                name_and_var( model_settings.dim ), []( auto x ) { return x == 1; }, basic_deff_msg );
         }
     }
 }
@@ -237,11 +266,8 @@ void print_settings( const SimulationOptions & options )
     fmt::print( "[Model]\n" );
     fmt::print( "    type {}\n", options.model_string );
 
-    // @TODO: Optionally print *all* settings to the console, including defaults that were set
-    if( options.model == Model::ActivityDrivenModel )
+    auto print_activity_settings = [&]( auto model_settings )
     {
-        auto model_settings = std::get<ActivityDrivenSettings>( options.model_settings );
-
         fmt::print( "    max_iterations {}\n", model_settings.max_iterations );
         fmt::print( "    dt {} \n", model_settings.dt );
         fmt::print( "    m {} \n", model_settings.m );
@@ -270,6 +296,20 @@ void print_settings( const SimulationOptions & options )
             fmt::print( "    reluctance_eps {}\n", model_settings.reluctance_eps );
             fmt::print( "    covariance_factor {}\n", model_settings.covariance_factor );
         }
+    };
+
+    // @TODO: Optionally print *all* settings to the console, including defaults that were set
+    if( options.model == Model::ActivityDrivenModel )
+    {
+        auto model_settings = std::get<ActivityDrivenSettings>( options.model_settings );
+
+        print_activity_settings( model_settings );
+    }
+    else if( options.model == Model::ActivityDrivenInertial )
+    {
+        auto model_settings = std::get<ActivityDrivenInertialSettings>( options.model_settings );
+        print_activity_settings( model_settings );
+        fmt::print( "    friction_coefficient {}\n", model_settings.friction_coefficient );
     }
     else if( options.model == Model::DeGroot )
     {
diff --git a/src/main.cpp b/src/main.cpp
index 9343730..3f429fa 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,6 +1,7 @@
 #include "config_parser.hpp"
 #include "models/DeGroot.hpp"
 #include "models/DeffuantModel.hpp"
+#include "models/InertialModel.hpp"
 #include "simulation.hpp"
 #include <fmt/format.h>
 #include <fmt/ostream.h>
@@ -70,6 +71,11 @@ int main( int argc, char * argv[] )
         simulation = std::make_unique<Seldon::Simulation<Seldon::ActivityDrivenModel::AgentT>>(
             simulation_options, network_file, agent_file );
     }
+    else if( simulation_options.model == Seldon::Config::Model::ActivityDrivenInertial )
+    {
+        simulation = std::make_unique<Seldon::Simulation<Seldon::InertialModel::AgentT>>(
+            simulation_options, network_file, agent_file );
+    }
     else if( simulation_options.model == Seldon::Config::Model::DeffuantModel )
     {
         simulation = std::make_unique<Seldon::Simulation<Seldon::DeffuantModel::AgentT>>(
diff --git a/src/models/ActivityDrivenModel.cpp b/src/models/ActivityDrivenModel.cpp
index 0a0ba06..0b0dd9c 100644
--- a/src/models/ActivityDrivenModel.cpp
+++ b/src/models/ActivityDrivenModel.cpp
@@ -1,6 +1,5 @@
 #include "models/ActivityDrivenModel.hpp"
 #include "network.hpp"
-#include "network_generation.hpp"
 #include "util/math.hpp"
 #include <cstddef>
 #include <random>
@@ -9,252 +8,8 @@
 namespace Seldon
 {
 
-ActivityDrivenModel::ActivityDrivenModel(
-    const Config::ActivityDrivenSettings & settings, NetworkT & network, std::mt19937 & gen )
-        : Model<ActivityDrivenModel::AgentT>( settings.max_iterations ),
-          network( network ),
-          contact_prob_list( std::vector<std::vector<NetworkT::WeightT>>( network.n_agents() ) ),
-          gen( gen ),
-          dt( settings.dt ),
-          m( settings.m ),
-          eps( settings.eps ),
-          gamma( settings.gamma ),
-          alpha( settings.alpha ),
-          homophily( settings.homophily ),
-          reciprocity( settings.reciprocity ),
-          K( settings.K ),
-          mean_activities( settings.mean_activities ),
-          mean_weights( settings.mean_weights ),
-          use_reluctances( settings.use_reluctances ),
-          reluctance_mean( settings.reluctance_mean ),
-          reluctance_sigma( settings.reluctance_sigma ),
-          reluctance_eps( settings.reluctance_eps ),
-          n_bots( settings.n_bots ),
-          bot_m( settings.bot_m ),
-          bot_activity( settings.bot_activity ),
-          bot_opinion( settings.bot_opinion ),
-          bot_homophily( settings.bot_homophily )
-{
-    get_agents_from_power_law();
-
-    if( mean_weights )
-    {
-        auto agents_copy = network.agents;
-        network          = NetworkGeneration::generate_fully_connected<AgentT>( network.n_agents() );
-        network.agents   = agents_copy;
-    }
-}
-
-double ActivityDrivenModel::homophily_weight( size_t idx_contacter, size_t idx_contacted )
-{
-    double homophily = this->homophily;
-
-    if( idx_contacted == idx_contacter )
-        return 0.0;
-
-    if( bot_present() && idx_contacter < n_bots )
-        homophily = this->bot_homophily[idx_contacter];
-
-    constexpr double tolerance = 1e-10;
-    auto opinion_diff
-        = std::abs( network.agents[idx_contacter].data.opinion - network.agents[idx_contacted].data.opinion );
-    opinion_diff = std::max( tolerance, opinion_diff );
-
-    return std::pow( opinion_diff, -homophily );
-}
-
-void ActivityDrivenModel::get_agents_from_power_law()
-{
-    std::uniform_real_distribution<> dis_opinion( -1, 1 ); // Opinion initial values
-    power_law_distribution<> dist_activity( eps, gamma );
-    truncated_normal_distribution<> dist_reluctance( reluctance_mean, reluctance_sigma, reluctance_eps );
-
-    auto mean_activity = dist_activity.mean();
-
-    // Initial conditions for the opinions, initialize to [-1,1]
-    // The activities should be drawn from a power law distribution
-    for( size_t i = 0; i < network.agents.size(); i++ )
-    {
-        network.agents[i].data.opinion = dis_opinion( gen ); // Draw the opinion value
-
-        if( !mean_activities )
-        { // Draw from a power law distribution (1-gamma)/(1-eps^(1-gamma)) * a^(-gamma)
-            network.agents[i].data.activity = dist_activity( gen );
-        }
-        else
-        {
-            network.agents[i].data.activity = mean_activity;
-        }
-
-        if( use_reluctances )
-        {
-            network.agents[i].data.reluctance = dist_reluctance( gen );
-        }
-    }
-
-    if( bot_present() )
-    {
-        for( size_t bot_idx = 0; bot_idx < n_bots; bot_idx++ )
-        {
-            network.agents[bot_idx].data.opinion  = bot_opinion[bot_idx];
-            network.agents[bot_idx].data.activity = bot_activity[bot_idx];
-        }
-    }
-}
-
-void ActivityDrivenModel::update_network_probabilistic()
-{
-    network.switch_direction_flag();
-
-    std::uniform_real_distribution<> dis_activation( 0.0, 1.0 );
-    std::uniform_real_distribution<> dis_reciprocation( 0.0, 1.0 );
-    std::vector<size_t> contacted_agents{};
-    reciprocal_edge_buffer.clear(); // Clear the reciprocal edge buffer
-    for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
-    {
-        // Test if the agent is activated
-        bool activated = dis_activation( gen ) < network.agents[idx_agent].data.activity;
-
-        if( activated )
-        {
-            // Implement the weight for the probability of agent `idx_agent` contacting agent `j`
-            // Not normalised since this is taken care of by the reservoir sampling
-
-            int m_temp = this->m;
-
-            if( bot_present() && idx_agent < n_bots )
-            {
-                m_temp = bot_m[idx_agent];
-            }
-
-            reservoir_sampling_A_ExpJ(
-                m_temp, network.n_agents(), [&]( int j ) { return homophily_weight( idx_agent, j ); }, contacted_agents,
-                gen );
-
-            // Fill the outgoing edges into the reciprocal edge buffer
-            for( const auto & idx_outgoing : contacted_agents )
-            {
-                reciprocal_edge_buffer.insert(
-                    { idx_agent, idx_outgoing } ); // insert the edge idx_agent -> idx_outgoing
-            }
-
-            // Set the *outgoing* edges
-            network.set_neighbours_and_weights( idx_agent, contacted_agents, 1.0 );
-        }
-        else
-        {
-            network.set_neighbours_and_weights( idx_agent, {}, {} );
-        }
-    }
-
-    // Reciprocity check
-    for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
-    {
-        // Get the outgoing edges
-        auto contacted_agents = network.get_neighbours( idx_agent );
-        // For each outgoing edge we check if the reverse edge already exists
-        for( const auto & idx_outgoing : contacted_agents )
-        {
-            // If the edge is not reciprocated
-            if( !reciprocal_edge_buffer.contains( { idx_outgoing, idx_agent } ) )
-            {
-                if( dis_reciprocation( gen ) < reciprocity )
-                {
-                    network.push_back_neighbour_and_weight( idx_outgoing, idx_agent, 1.0 );
-                }
-            }
-        }
-    }
-
-    network.toggle_incoming_outgoing(); // switch direction, so that we have incoming edges
-}
-
-void ActivityDrivenModel::update_network_mean()
-{
-    using WeightT = NetworkT::WeightT;
-    std::vector<WeightT> weights( network.n_agents(), 0.0 );
-
-    // Set all weights to zero in the beginning
-    for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
-    {
-        network.set_weights( idx_agent, weights );
-        contact_prob_list[idx_agent] = weights; // set to zero
-    }
-
-    auto probability_helper = []( double omega, size_t m )
-    {
-        double p = 0;
-        for( size_t i = 1; i <= m; i++ )
-            p += ( std::pow( -omega, i + 1 ) + omega ) / ( omega + 1 );
-        return p;
-    };
-
-    for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
-    {
-        // Implement the weight for the probability of agent `idx_agent` contacting agent `j`
-        // Not normalised since this is taken care of by the reservoir sampling
-
-        double normalization = 0;
-        for( size_t k = 0; k < network.n_agents(); k++ )
-        {
-            normalization += homophily_weight( idx_agent, k );
-        }
-
-        // Go through all the neighbours of idx_agent
-        // Calculate the probability of i contacting j (in 1 to m rounds, assuming
-        // the agent is activated
-        int m_temp = m;
-        if( bot_present() && idx_agent < n_bots )
-        {
-            m_temp = bot_m[idx_agent];
-        }
-
-        double activity = std::max( 1.0, network.agents[idx_agent].data.activity );
-        for( size_t j = 0; j < network.n_agents(); j++ )
-        {
-            double omega = homophily_weight( idx_agent, j ) / normalization;
-            // We can calculate the probability of i contacting j ( i->j )
-            // Update contact prob_list (outgoing)
-            contact_prob_list[idx_agent][j] = activity * probability_helper( omega, m_temp );
-        }
-    }
-
-    for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
-    {
-        // Calculate the actual weights and reciprocity
-        for( size_t j = 0; j < network.n_agents(); j++ )
-        {
-            double prob_contact_ij = contact_prob_list[idx_agent][j]; // outgoing probabilites
-            double prob_contact_ji = contact_prob_list[j][idx_agent];
-
-            // Set the incoming agent weight, j-i in weight list
-            double & win_ji = network.get_weights( j )[idx_agent];
-            win_ji += prob_contact_ij;
-
-            // Handle the reciprocity for j->i
-            // Update incoming weight i-j
-            double & win_ij = network.get_weights( idx_agent )[j];
-
-            // The probability of reciprocating is
-            win_ij += ( 1.0 - prob_contact_ji ) * reciprocity * prob_contact_ij;
-        }
-    }
-}
-
-void ActivityDrivenModel::update_network()
-{
-
-    if( !mean_weights )
-    {
-        update_network_probabilistic();
-    }
-    else
-    {
-        update_network_mean();
-    }
-}
-
-void ActivityDrivenModel::iteration()
+template<>
+void ActivityDrivenModelAbstract<ActivityAgent>::iteration()
 {
     Model<AgentT>::iteration();
 
@@ -265,19 +20,21 @@ void ActivityDrivenModel::iteration()
     get_euler_slopes( k1_buffer, [this]( size_t i ) { return network.agents[i].data.opinion; } );
     // k_2  =   hf(x_n+1/2h,y_n+1/2k_1)
     get_euler_slopes(
-        k2_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + 0.5 * this->k1_buffer[i]; } );
+        k2_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + 0.5 * dt * this->k1_buffer[i]; } );
     // k_3  =   hf(x_n+1/2h,y_n+1/2k_2)
     get_euler_slopes(
-        k3_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + 0.5 * this->k2_buffer[i]; } );
+        k3_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + 0.5 * dt * this->k2_buffer[i]; } );
     // k_4  =   hf(x_n+h,y_n+k_3)
-    get_euler_slopes( k4_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + this->k3_buffer[i]; } );
+    get_euler_slopes(
+        k4_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + dt * this->k3_buffer[i]; } );
 
     // Update the agent opinions
     for( size_t idx_agent = 0; idx_agent < network.n_agents(); ++idx_agent )
     {
         // y_(n+1) =   y_n+1/6k_1+1/3k_2+1/3k_3+1/6k_4+O(h^5)
         network.agents[idx_agent].data.opinion
-            += ( k1_buffer[idx_agent] + 2 * k2_buffer[idx_agent] + 2 * k3_buffer[idx_agent] + k4_buffer[idx_agent] )
+            += dt
+               * ( k1_buffer[idx_agent] + 2 * k2_buffer[idx_agent] + 2 * k3_buffer[idx_agent] + k4_buffer[idx_agent] )
                / 6.0;
     }
 
@@ -289,4 +46,6 @@ void ActivityDrivenModel::iteration()
         }
     }
 }
+
+template class ActivityDrivenModelAbstract<ActivityAgent>;
 } // namespace Seldon
\ No newline at end of file
diff --git a/src/models/InertialModel.cpp b/src/models/InertialModel.cpp
new file mode 100644
index 0000000..df46cc1
--- /dev/null
+++ b/src/models/InertialModel.cpp
@@ -0,0 +1,77 @@
+#include "models/InertialModel.hpp"
+#include "agents/inertial_agent.hpp"
+#include "network.hpp"
+#include "util/math.hpp"
+#include <cstddef>
+#include <random>
+#include <vector>
+
+namespace Seldon
+{
+
+InertialModel::InertialModel(
+    const Config::ActivityDrivenInertialSettings & settings, NetworkT & network, std::mt19937 & gen )
+        : ActivityDrivenModelAbstract<InertialAgent>( settings, network, gen ),
+          friction_coefficient( settings.friction_coefficient )
+{
+}
+
+// X(t+dt)
+// Also updates the position
+void InertialModel::calc_position()
+{
+    // Calculating 'drift' = a(t)-friction
+    get_euler_slopes( drift_t_buffer, [this]( size_t i ) { return network.agents[i].data.opinion; } );
+
+    for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
+    {
+        auto & agent_data    = network.agents[idx_agent].data;
+        auto accleration     = drift_t_buffer[idx_agent] - friction_coefficient * agent_data.velocity;
+        double next_position = agent_data.opinion + agent_data.velocity * dt + 0.5 * (accleration)*dt * dt;
+
+        // Update the position to the new position
+        agent_data.opinion = next_position;
+    }
+}
+
+// V(t+dt)
+// Also updates the velocity
+void InertialModel::calc_velocity()
+{
+    // Calculating new 'drift'
+    get_euler_slopes( drift_next_t_buffer, [this]( size_t i ) { return network.agents[i].data.opinion; } );
+
+    for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
+    {
+        auto & agent_data    = network.agents[idx_agent].data;
+        double next_velocity = agent_data.velocity
+                               + 0.5 * dt
+                                     * ( drift_t_buffer[idx_agent] - friction_coefficient * agent_data.velocity
+                                         + drift_next_t_buffer[idx_agent] );
+        next_velocity /= 1.0 + 0.5 * friction_coefficient * dt;
+
+        // Update velocity
+        agent_data.velocity = next_velocity;
+    }
+}
+
+void InertialModel::iteration()
+{
+    Model<AgentT>::iteration();
+
+    update_network();
+
+    // Use Velocity Verlet algorithm
+    calc_position();
+    calc_velocity();
+
+    if( bot_present() )
+    {
+        for( size_t bot_idx = 0; bot_idx < n_bots; bot_idx++ )
+        {
+            network.agents[bot_idx].data.opinion = bot_opinion[bot_idx];
+        }
+    }
+}
+
+} // namespace Seldon
\ No newline at end of file
diff --git a/test/res/1bot_1agent_inertial.toml b/test/res/1bot_1agent_inertial.toml
new file mode 100644
index 0000000..6520309
--- /dev/null
+++ b/test/res/1bot_1agent_inertial.toml
@@ -0,0 +1,42 @@
+[simulation]
+model = "ActivityDrivenInertial"
+rng_seed = 120           # Leaving this empty will pick a random seed
+
+[io]
+# n_output_network = 1  # Write the network every 20 iterations
+# n_output_agents = 1
+print_progress = false # Print the iteration time ; if not set, then does not print
+
+[model]
+max_iterations = 1000 # If not set, max iterations is infinite
+
+[ActivityDrivenInertial]
+dt = 0.001              # Timestep for the integration of the coupled ODEs
+m = 1                   # Number of agents contacted, when the agent is active
+eps = 1                 # Minimum activity epsilon; a_i belongs to [epsilon,1]
+gamma = 2.1             # Exponent of activity power law distribution of activities
+reciprocity = 1         # probability that when agent i contacts j via weighted reservoir sampling, j also sends feedback to i. So every agent can have more than m incoming connections
+homophily = 0.5         # aka beta. if zero, agents pick their interaction partners at random
+alpha = 1.5             # Controversialness of the issue, must be greater than 0.
+K = 2.0                 # Social interaction strength
+mean_activities = false # Use the mean value of the powerlaw distribution for the activities of all agents
+mean_weights = false    # Use the meanfield approximation of the network edges
+
+reluctances = true     # Assigns a "reluctance" (m_i) to each agent. By default; false and every agent has a reluctance of 1
+reluctance_mean = 1.5  # Mean of distribution before drawing from a truncated normal distribution (default set to 1.0)
+reluctance_sigma = 0.1 # Width of normal distribution (before truncating)
+reluctance_eps = 0.01  # Minimum such that the normal distribution is truncated at this value
+
+n_bots = 1 # The number of bots to be used; if not specified defaults to 0 (which means bots are deactivated)
+# Bots are agents with fixed opinions and different parameters, the parameters are specified in the following lists
+# If n_bots is smaller than the length of any of the lists, the first n_bots entries are used. If n_bots is greater the code will throw an exception.
+bot_m = [1]           # If not specified, defaults to `m`
+bot_homophily = [0.7] # If not specified, defaults to `homophily`
+bot_activity = [1.0]  # If not specified, defaults to 0
+bot_opinion = [2]     # The fixed opinions of the bots
+
+friction_coefficient = 0.5 # Friction coefficient 
+
+[network]
+number_of_agents = 2
+connections_per_agent = 1
diff --git a/test/test_activity_inertial.cpp b/test/test_activity_inertial.cpp
new file mode 100644
index 0000000..60f3eef
--- /dev/null
+++ b/test/test_activity_inertial.cpp
@@ -0,0 +1,76 @@
+#include "catch2/matchers/catch_matchers.hpp"
+#include "models/InertialModel.hpp"
+#include "util/math.hpp"
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <complex>
+#include <config_parser.hpp>
+#include <filesystem>
+#include <network_io.hpp>
+#include <optional>
+#include <simulation.hpp>
+
+namespace fs = std::filesystem;
+
+TEST_CASE( "Test the probabilistic inertial activity driven model with one bot and one agent", "[inertial1Bot1Agent]" )
+{
+    using namespace Seldon;
+    using namespace Catch::Matchers;
+    using AgentT = InertialModel::AgentT;
+
+    auto proj_root_path = fs::current_path();
+    auto input_file     = proj_root_path / fs::path( "test/res/1bot_1agent_inertial.toml" );
+
+    auto options = Config::parse_config_file( input_file.string() );
+    Config::print_settings( options );
+
+    auto simulation = Simulation<AgentT>( options, std::nullopt, std::nullopt );
+
+    // We need an output path for Simulation, but we won't write anything out there
+    fs::path output_dir_path = proj_root_path / fs::path( "test/output_inertial" );
+
+    // Get the bot opinion (which won't change)
+    auto bot   = simulation.network.agents[0];
+    auto x_bot = bot.data.opinion; // Bot opinion
+    fmt::print( "We have set the bot opinion = {}\n", x_bot );
+
+    // Get the initial agent opinion
+    auto & agent = simulation.network.agents[1];
+    auto x_0     = agent.data.opinion; // Agent opinion
+    fmt::print( "We have set agent x_0 = {}\n", x_0 );
+
+    simulation.run( output_dir_path );
+
+    auto model_settings     = std::get<Seldon::Config::ActivityDrivenInertialSettings>( options.model_settings );
+    auto K                  = model_settings.K;
+    auto alpha              = model_settings.alpha;
+    auto iterations         = model_settings.max_iterations.value();
+    auto dt                 = model_settings.dt;
+    std::complex<double> mu = model_settings.friction_coefficient;
+    auto time_elapsed       = iterations * dt;
+
+    // Final agent and bot opinions after the simulation run
+    auto x_t        = agent.data.opinion;
+    auto x_t_bot    = bot.data.opinion;
+    auto reluctance = agent.data.reluctance;
+
+    fmt::print( "Agent reluctance is = {}\n", reluctance );
+
+    // The bot opinion should not change during the simulation
+    REQUIRE_THAT( x_t_bot, WithinAbs( x_bot, 1e-16 ) );
+
+    // C = K/m tanh (alpha*x_bot)
+    auto C = K / reluctance * std::tanh( alpha * x_bot );
+
+    fmt::print( "C =     {}\n", C );
+    auto a1 = 0.5 * ( -std::sqrt( mu * mu - 4.0 ) - mu );
+    auto a2 = 0.5 * ( std::sqrt( mu * mu - 4.0 ) - mu );
+    auto c1 = ( x_0 - C ) / ( 1.0 - a1 / a2 );
+    auto c2 = -c1 * a1 / a2;
+    // Test that the agent opinion matches the analytical solution for an agent with a bot
+    // Analytical solution is
+    // x_t = c1 * exp(a1*t) + c2 *exp(a2*t) + C
+    auto x_t_analytical = c1 * std::exp( a1 * time_elapsed ) + c2 * std::exp( a2 * time_elapsed ) + C;
+
+    REQUIRE_THAT( x_t, WithinAbs( x_t_analytical.real(), 1e-5 ) );
+}
\ No newline at end of file
diff --git a/test/test_network.cpp b/test/test_network.cpp
index f1d2bfa..57239e7 100644
--- a/test/test_network.cpp
+++ b/test/test_network.cpp
@@ -112,7 +112,7 @@ TEST_CASE( "Testing the network class" )
                 auto weight    = buffer_w[i_neighbour];
                 std::tuple<size_t, size_t, Network::WeightT> edge{
                     neighbour, i_agent, weight
-                }; // Note that i_agent and neighbour are flipped compared to before
+                };                                     // Note that i_agent and neighbour are flipped compared to before
                 REQUIRE( old_edges.contains( edge ) ); // can we find the transposed edge?
                 old_edges.extract( edge );             // extract the edge afterwards
             }