diff --git a/examples/quickstart-fr.py b/examples/quickstart-fr.py
new file mode 100644
index 00000000..5d0d84ba
--- /dev/null
+++ b/examples/quickstart-fr.py
@@ -0,0 +1,78 @@
+import os
+import dotenv
+import logging
+
+import mobility
+
+import polars as pl
+
+dotenv.load_dotenv()
+
+
+mobility.set_params(
+    # package_data_folder_path=os.environ["MOBILITY_PACKAGE_DATA_FOLDER"],
+    # project_data_folder_path=os.environ["MOBILITY_PROJECT_DATA_FOLDER"]
+        package_data_folder_path="D:/mobility-data",
+        project_data_folder_path="D:/test-09",
+        debug=False
+)
+
+# Using Foix (a small town) and a limited radius for quick results
+transport_zones = mobility.TransportZones("fr-09122", radius = 10)
+
+# Using EMP, the latest national mobility survey for France
+emp = mobility.EMPMobilitySurvey()
+
+# Creating a synthetic population of 1000 for the area
+pop = mobility.Population(transport_zones, sample_size = 1000)
+
+modes = [mobility.CarMode(transport_zones), mobility.WalkMode(transport_zones), mobility.BicycleMode(transport_zones)]
+surveys = [emp]
+motives = [mobility.HomeMotive(), mobility.WorkMotive(), mobility.OtherMotive(population=pop)]
+
+
+# Simulating the trips for this population for three modes : car, walk and bicyle, and only home and work motives (OtherMotive is mandatory)
+pop_trips = mobility.PopulationTrips(
+    pop,
+    modes,
+    motives,
+    surveys,
+    parameters=mobility.PopulationTripsParameters(n_iterations=4, k_mode_sequences=42)
+    )
+
+pop_trips_base = mobility.PopulationTrips(
+    pop,
+    modes,
+    motives,
+    surveys,
+    parameters=mobility.PopulationTripsParameters(n_iterations=4, seed=49283092, k_mode_sequences=42)
+    )
+
+tz = transport_zones.get()
+
+
+# You can get the weekday trips to inspect them
+weekday_flows = pop_trips.get()["weekday_flows"].collect()
+weekday_flows_base = pop_trips_base.get()["weekday_flows"].collect()
+
+# You can also plot the flows, with labels for the cities that are bigger than their neighbours
+labels=pop_trips.get_prominent_cities()
+
+# You can can also get global metrics that will be compared to the theoetical values for this population 
+cost_per_person = pop_trips.evaluate("cost_per_person", plot_delta=True, compare_with=pop_trips_base, labels=labels, plot=True)
+dist_per_person = pop_trips.evaluate("distance_per_person", plot_delta=True, compare_with=pop_trips_base, plot=True)
+time_per_person = pop_trips.evaluate("time_per_person", plot_delta=True, compare_with=pop_trips_base, plot=True)
+ghg_per_person = pop_trips.evaluate("ghg_per_person", plot_delta=True, compare_with=pop_trips_base, plot=True)
+global_metrics = pop_trips.evaluate("global_metrics")
+
+
+pop_trips.plot_modal_share(mode="public_transport", labels=labels)
+pop_trips.plot_modal_share(mode="bicycle")
+pop_trips.plot_modal_share(mode="walk")
+cms = pop_trips.plot_modal_share(mode="car")
+
+# OD flows between transport zones
+pop_trips.plot_od_flows(mode="car", level_of_detail=1, labels=labels)
+pop_trips.plot_od_flows(mode="walk", level_of_detail=1)
+pop_trips.plot_od_flows(mode="bicycle", level_of_detail=1)
+pop_trips.plot_od_flows(mode="public_transport")
\ No newline at end of file
diff --git a/examples/quickstart-fr.py.py b/examples/quickstart-fr.py.py
deleted file mode 100644
index ebc32e82..00000000
--- a/examples/quickstart-fr.py.py
+++ /dev/null
@@ -1,40 +0,0 @@
-import os
-import dotenv
-
-import mobility
-
-dotenv.load_dotenv()
-
-mobility.set_params(
-    # package_data_folder_path=os.environ["MOBILITY_PACKAGE_DATA_FOLDER"],
-    # project_data_folder_path=os.environ["MOBILITY_PROJECT_DATA_FOLDER"]
-        package_data_folder_path="D:/mobility-data",
-        project_data_folder_path="D:/test-09",
-)
-
-# Using Foix (a small town) and a limited radius for quick results
-transport_zones = mobility.TransportZones("fr-09122", radius = 10)
-
-# Using EMP, the latest national mobility survey for France
-emp = mobility.EMPMobilitySurvey()
-
-# Creating a synthetic population of 1000 for the area
-pop = mobility.Population(transport_zones, sample_size = 1000)
-
-# Simulating the trips for this population for three modes : car, walk and bicyle, and only home and work motives (OtherMotive is mandatory)
-pop_trips = mobility.PopulationTrips(
-    pop,
-    [mobility.CarMode(transport_zones), mobility.WalkMode(transport_zones), mobility.BicycleMode(transport_zones)],
-    [mobility.HomeMotive(), mobility.WorkMotive(), mobility.OtherMotive(population=pop)],
-    [emp]
-    )
-
-# You can get the weekday trips to inspect them
-weekday_flows = pop_trips.get()["weekday_flows"].collect()
-
-# You can can also get global metrics that will be compared to the theoetical values for this population 
-global_metrics = pop_trips.evaluate("global_metrics")
-
-# You can also plot the flows, with labels for the cities that are bigger than their neighbours
-labels = pop_trips.get_prominent_cities()
-pop_trips.plot_od_flows(labels=labels)
diff --git a/mobility/__init__.py b/mobility/__init__.py
index 8a446432..12126389 100644
--- a/mobility/__init__.py
+++ b/mobility/__init__.py
@@ -1,5 +1,7 @@
 from .set_params import set_params
 
+from .model_parameters import Parameter, ParameterSet
+
 from .study_area import StudyArea
 from .transport_zones import TransportZones
 
@@ -31,7 +33,7 @@
 
 from .cost_of_time_parameters import CostOfTimeParameters
 
-from .choice_models.transport_mode_choice_model import TransportModeChoiceModel
+#from .choice_models.transport_mode_choice_model import TransportModeChoiceModel
 from .parsers import LocalAdminUnits
 
 
diff --git a/mobility/choice_models/destination_sequence_sampler.py b/mobility/choice_models/destination_sequence_sampler.py
index 13506264..5246032c 100644
--- a/mobility/choice_models/destination_sequence_sampler.py
+++ b/mobility/choice_models/destination_sequence_sampler.py
@@ -57,13 +57,13 @@ def run(
             transport_zones,
             remaining_sinks,
             costs,
-            parameters.cost_uncertainty_sd
+            parameters["cost_uncertainty_sd"]
         )
         
         dest_prob = self.get_destination_probability(
             utilities,
             motives,
-            parameters.dest_prob_cutoff
+            parameters["dest_prob_cutoff"]
         )
 
         chains = (
@@ -74,7 +74,7 @@ def run(
         )
         
         chains = self.spatialize_anchor_motives(chains, dest_prob, seed)
-        chains = self.spatialize_other_motives(chains, dest_prob, costs, parameters.alpha, seed)
+        chains = self.spatialize_other_motives(chains, dest_prob, costs, parameters["alpha"], seed)
         
         dest_sequences = ( 
             chains
diff --git a/mobility/choice_models/population_trips.py b/mobility/choice_models/population_trips.py
index 8b916d83..80a0ebd3 100644
--- a/mobility/choice_models/population_trips.py
+++ b/mobility/choice_models/population_trips.py
@@ -151,7 +151,7 @@ def __init__(
             "motives": motives,
             "modes": modes,
             "surveys": surveys,
-            "parameters": parameters
+            "parameters": parameters.to_dict()
         }
         
         project_folder = pathlib.Path(os.environ["MOBILITY_PROJECT_DATA_FOLDER"])
@@ -173,7 +173,9 @@ def __init__(
         }
         
         super().__init__(inputs, cache_path)
-        
+    
+    def has_multiple_seeds(self):
+        return (len(self.parameters.seeds) > 1)
         
     def resolve_parameters(
             self,
@@ -299,7 +301,7 @@ def create_and_get_asset(self):
         
         demand_groups.write_parquet(self.cache_path["demand_groups"])
         
-        if self.parameters.simulate_weekend:
+        if self.parameters["simulate_weekend"]:
             
             weekend_flows, weekend_sinks, demand_groups, weekend_costs, weekend_chains = self.run_model(is_weekday=False)
         
@@ -370,7 +372,7 @@ def run_model(self, is_weekday):
         
         remaining_sinks = sinks.clone()
         
-        for iteration in range(1, parameters.n_iterations+1):
+        for iteration in range(1, parameters["n_iterations"]+1):
             
             logging.info(f"Iteration n°{iteration}")
             
@@ -421,7 +423,7 @@ def run_model(self, is_weekday):
             costs = self.state_updater.get_new_costs(
                 costs,
                 iteration,
-                parameters.n_iter_per_cost_update,
+                parameters["n_iter_per_cost_update"],
                 current_states_steps,
                 costs_aggregator
             )
@@ -434,7 +436,7 @@ def run_model(self, is_weekday):
             
             
         costs = costs_aggregator.get_costs_by_od_and_mode(
-            ["cost", "distance", "time", "ghg_emissions"],
+            ["distance", "time", "ghg_emissions", "cost"],
             congestion=True
         )
         
@@ -539,7 +541,7 @@ def evaluate(self, metric, **kwargs):
         return evaluation
         
 
-    def plot_modal_share(self, zone="origin", mode="car", period="weekdays",
+    def plot_modal_share(self, zone="origin", mode="car", period="weekdays", plot=True,
                          labels=None, labels_size=[10, 6, 4], labels_color="black"):
         """
         Plot modal share for the given mode in the origin or destination zones during weekdays or weekends.
@@ -560,7 +562,6 @@ def plot_modal_share(self, zone="origin", mode="car", period="weekdays",
             Mode share for the given mode in each transport zone.
 
         """
-        logging.info(f"🗺️ Plotting {mode} modal share for {zone} zones during {period}")
 
         if period == "weekdays":
             population_df = self.get()["weekday_flows"].collect().to_pandas()
@@ -588,24 +589,27 @@ def plot_modal_share(self, zone="origin", mode="car", period="weekdays",
             mode_name = mode.capitalize()
         mode_share = mode_share[mode_share["mode"] == mode]
 
-        transport_zones_df = self.population.transport_zones.get()
-        
-        gc = gpd.GeoDataFrame(
-            transport_zones_df.merge(mode_share, how="left", right_on=left_column, left_on="transport_zone_id", suffixes=('', '_z'))).fillna(0)
-        gcp = gc.plot("modal_share", legend=True)
-        gcp.set_axis_off()
-        plt.title(f"{mode_name} share per {zone} transport zone ({period})")
+        if plot:
+            logging.info(f"🗺️ Plotting {mode} modal share for {zone} zones during {period}")
 
-        if isinstance(labels, gpd.GeoDataFrame):
-            self.__show_labels__(labels, labels_size, labels_color)        
-        
-        plt.show()
+            transport_zones_df = self.population.transport_zones.get()
+            
+            gc = gpd.GeoDataFrame(
+                transport_zones_df.merge(mode_share, how="left", right_on=left_column, left_on="transport_zone_id", suffixes=('', '_z'))).fillna(0)
+            gcp = gc.plot("modal_share", legend=True)
+            gcp.set_axis_off()
+            plt.title(f"{mode_name} share per {zone} transport zone ({period})")
+    
+            if isinstance(labels, gpd.GeoDataFrame):
+                self.__show_labels__(labels, labels_size, labels_color)        
+            
+            plt.show()
 
         return mode_share
 
     def plot_od_flows(self, mode="all", motive="all", period="weekdays", level_of_detail=1,
                       n_largest=2000, color="blue", transparency=0.2, zones_color="xkcd:light grey",
-                      labels=None, labels_size=[10, 6, 4], labels_color="black"):
+                      labels=None, labels_size=[10, 6, 4], labels_color="black", save=False):
         """
         Plot flows between the different zones for the given mode, motive, period and level of detail.
 
@@ -685,8 +689,8 @@ def plot_od_flows(self, mode="all", motive="all", period="weekdays", level_of_de
         # Put a legend for width on bottom right, title on the top
         x_min = float(biggest_flows[["x"]].min().iloc[0])
         y_min = float(biggest_flows[["y"]].min().iloc[0])
-        plt.plot([x_min, x_min+4000], [y_min, y_min], linewidth=2, color=color)
-        plt.text(x_min+6000, y_min-1000, "1 000", color=color)
+        plt.plot([x_min-10000, x_min-8000], [y_min-2000, y_min-2000], linewidth=2, color=color)
+        plt.text(x_min-6000, y_min-3000, "1 000", color=color)
         plt.title(f"{mode_name} flows between transport zones on {period}")
 
         # Draw all origin-destinations
@@ -696,6 +700,10 @@ def plot_od_flows(self, mode="all", motive="all", period="weekdays", level_of_de
 
         if isinstance(labels, gpd.GeoDataFrame):
             self.__show_labels__(labels, labels_size, labels_color)
+            
+        if save:
+            text = "plot-" + str(self.parameters["seed"]) + ".png"
+            plt.savefig(text)
 
         plt.show()
 
@@ -757,13 +765,14 @@ def get_prominent_cities(self, n_cities=20, n_levels=3, distance_km=2):
 
         # If an admin unit is too close to a bigger admin unit, make it less prominent
         for i in range(n_cities//2):
-            coords = flows_per_commune.loc[i, "geometry"]
-            geoflows["dists"] = geoflows["geometry"].distance(coords)
-            # Use distance_km here
-            geoflows.loc[
-                ((geoflows["dists"] < distance_km*1000) & (geoflows.index > i)), "prominence"
-                ] = geoflows["prominence"] + 2
-            geoflows = geoflows.sort_values(by="prominence").reset_index(drop=True)
+            if i < len(flows_per_commune):
+                coords = flows_per_commune.loc[i, "geometry"]
+                geoflows["dists"] = geoflows["geometry"].distance(coords)
+                # Use distance_km here
+                geoflows.loc[
+                    ((geoflows["dists"] < distance_km*1000) & (geoflows.index > i)), "prominence"
+                    ] = geoflows["prominence"] + 2
+                geoflows = geoflows.sort_values(by="prominence").reset_index(drop=True)
 
         # Keep only the most prominent admin units and add their centroids
         geoflows = geoflows[geoflows["prominence"] <= n_levels]
diff --git a/mobility/choice_models/population_trips_multi.py b/mobility/choice_models/population_trips_multi.py
new file mode 100644
index 00000000..c7f1a08b
--- /dev/null
+++ b/mobility/choice_models/population_trips_multi.py
@@ -0,0 +1,235 @@
+import os
+import pathlib
+import logging
+import shutil
+import random
+import warnings
+
+import geopandas as gpd
+import matplotlib.pyplot as plt
+import polars as pl
+
+import plotly.express as px
+import plotly.colors as pc
+
+
+from typing import List
+
+from mobility.population import Population
+from mobility.choice_models.population_trips import PopulationTrips
+from mobility.choice_models.population_trips_parameters import PopulationTripsParameters
+from mobility.choice_models.results import Results
+from mobility.motives import Motive
+from mobility.transport_modes.transport_mode import TransportMode
+from mobility.parsers.mobility_survey import MobilitySurvey
+
+#New library
+from PIL import Image
+
+# Script ones
+import dotenv
+import mobility
+
+class MultiSeedPopulationTrips:
+    """Class created to manage multi-seed PopulationTrips. Parameters are the same than for PopulationTrips. """
+    
+    def __init__(
+            
+            self,
+            population: Population,
+            modes: List[TransportMode] = None,
+            motives: List[Motive] = None,
+            surveys: List[MobilitySurvey] = None,
+            parameters: PopulationTripsParameters = None,
+            seeds: List[int] = [0]
+            
+        ):
+        
+        list_population_trips = []
+        
+        for seed in seeds:
+
+            print("init seed", seed)
+            parameters_seed = parameters
+            parameters_seed.seed = seed
+            pop_seed = PopulationTrips(
+                population = population,
+                modes = modes,
+                motives = motives,
+                surveys=surveys,
+                parameters = parameters_seed
+                )
+            list_population_trips.append(pop_seed)
+        
+        self.trips = list_population_trips
+        
+    def get(self):
+        return self.trips
+    
+    def evaluate(self, metric, plot_diff=False, plot=True, compare_with=None, variable=None, **kwargs):
+        if metric in ["distance_per_person", "time_per_person", "cost_per_person"]: #todo: , "ghg_per_person"
+            metric_per_person = metric
+            metric = metric.replace("_per_person", "")
+            
+            metric_all = self.trips[0].evaluate(metric_per_person).select("transport_zone_id", metric, "n_persons").group_by("transport_zone_id").agg(pl.col(metric).sum(), pl.col("n_persons").sum())
+            i = 1
+            for pop in self.trips[1:]:
+                metric_add = pop.evaluate(metric_per_person).select("transport_zone_id", metric, "n_persons").group_by("transport_zone_id").agg(pl.col(metric).sum(), pl.col("n_persons").sum())
+                suffix = "_" + str(i)
+                metric_all = metric_all.join(metric_add, on="transport_zone_id", suffix=suffix)
+                i += 1
+            metric_all = metric_all.sort(by=metric)
+            regex = "^" + metric + ".*$"
+            metric_all = metric_all.with_columns(metric_mean=pl.mean_horizontal(regex),
+                                                 n_persons_mean=pl.mean_horizontal("^n_persons.*$"))
+            metric_all = metric_all.with_columns(metric_per_person=pl.col("metric_mean")/pl.col("n_persons_mean"))
+            
+            if plot_diff:
+                logging.info(f"Plotting {metric}, each curve being a different seed. Curves should be close.")
+                pandas_metric=metric_all.select(pl.col(regex)).to_pandas()
+                pandas_metric.plot()
+                
+            if plot:
+                pass #todo : add the mean and median map for the metric
+                
+            if compare_with is not None:
+                metric_all = metric_all.join(compare_with.select("transport_zone_id", "metric_per_person"), on="transport_zone_id", suffix="_base")
+                metric_all = metric_all.with_columns(delta=pl.col("metric_per_person")-pl.col("metric_per_person_base"))
+
+            return metric_all
+        
+        elif metric == "metrics_by_variable":
+            print(variable)
+            metric_all = self.trips[0].evaluate("metrics_by_variable", variable=variable, plot=False)
+            i = 1
+            for pop in self.trips[1:]:
+                metric_add = pop.evaluate("metrics_by_variable", variable=variable, plot=False).select("variable", variable, "value", "delta", "delta_relative")
+                suffix = "_" + str(i)
+                metric_all = metric_all.join(metric_add, on=("variable", variable), suffix=suffix)
+                i += 1
+                
+            if plot:              
+                comparison_plot_df = (
+                    metric_all
+                    .rename({"value":"value_0"})
+                    .select(["variable", variable, pl.col("^value.*$")])
+                    .melt(["variable", variable], variable_name="value_type")
+                    .sort(variable)
+                )
+                
+                value_types = comparison_plot_df["value_type"].unique()
+                value_types = [v for v in value_types if v != "value_ref"] + ["value_ref"]
+
+                
+                blues = pc.sequential.Blues_r
+                blue_values = [v for v in value_types if v != "value_ref"]
+                
+                color_map = {
+                    vt: blues[i % len(blues)]
+                    for i, vt in enumerate(blue_values)
+                }
+                color_map["value_ref"] = "red"
+                
+                fig = px.bar(
+                    comparison_plot_df,
+                    x=variable,
+                    y="value",
+                    color="value_type",
+                    facet_col="variable",
+                    barmode="group",
+                    facet_col_spacing=0.05,
+                    color_discrete_map=color_map,
+                    category_orders={"value_type": value_types}
+                )
+                fig.update_yaxes(matches=None, showticklabels=True)
+                fig.show("browser")
+            
+        else:
+            return NotImplemented
+    
+            
+    def plot_modal_share(self, plot_mode="all", **kwargs):
+        if plot_mode == "all":
+            logging.info("Plotting all modal shares for the given seeds")
+            for t in self.trips:
+                t.plot_modal_share(**kwargs)
+        else:
+            return NotImplemented
+    
+    def plot_od_flows(self, plot_mode="all", **kwargs):
+        """
+        Plots flows for the given seeds. Can show them all or create a GIF.
+
+        Parameters
+        ----------
+        plot_mode : TYPE, optional
+            Plot mode : "all" shows OD flows for all seeds. "gif" creates a gif. The default is "all".
+        **kwargs : TYPE
+            Same parameters than the function plot_od_flows from PopulationTrips
+        """
+        if plot_mode == "all":
+            logging.info("Plotting all OD flows for the given seeds")
+            for t in self.trips:
+                t.plot_od_flows(**kwargs)
+        elif plot_mode == "gif":
+            frames = []
+            for t in self.trips:
+                t.plot_od_flows(save=True, **kwargs)
+                text = "plot-" + str(t.parameters["seed"]) + ".png"
+                frames.append(Image.open(text))
+            
+        frames[0].save(
+            "output.gif",
+            format="GIF",
+            append_images=frames[1:],
+            save_all=True,
+            duration=400,  # milliseconds per frame
+            loop=0         # 0 = infinite loop
+        )    
+        
+    def get_prominent_cities(self, **kwargs):
+        """Returning the cities of the first simulation as it is consistent over simulations"""
+        return self.trips[0].get_prominent_cities(**kwargs)
+        
+if __name__ == "__main__":
+    dotenv.load_dotenv()
+    os.environ["MOBILITY_GTFS_DOWNLOAD_DATE"] = "2025/12/01"
+    mobility.set_params(
+            package_data_folder_path="D:/mobility-data",
+            project_data_folder_path="D:/test-09",
+            debug=False
+    )
+
+    transport_zones = mobility.TransportZones("fr-29019", radius = 40, level_of_detail=1)
+    emp = mobility.EMPMobilitySurvey()
+    pop = mobility.Population(transport_zones, sample_size = 1000)
+    modes = [mobility.CarMode(transport_zones), mobility.WalkMode(transport_zones), mobility.BicycleMode(transport_zones), mobility.PublicTransportMode(transport_zones)]
+    surveys = [emp]
+    motives = [mobility.HomeMotive(), mobility.WorkMotive(), mobility.OtherMotive(population=pop)]
+
+    # Simulating the trips for this population for three modes : car, walk and bicyle, and only home and work motives (OtherMotive is mandatory)
+    pop_trips = MultiSeedPopulationTrips(
+        pop,
+        modes,
+        motives,
+        surveys,
+        parameters=mobility.PopulationTripsParameters(n_iterations=4, k_mode_sequences=6),
+        seeds = [0, 42, 69, 99, 123]
+        )
+        
+    # for t in pop_trips.get():
+    #     g = t.get()
+    
+    #labels = pop_trips.get_prominent_cities()
+    #pop_trips.plot_modal_share(mode="walk")
+    #pop_trips.plot_od_flows(mode="gif", labels=labels)
+    d = pop_trips.evaluate("distance_per_person")
+    t = pop_trips.evaluate("time_per_person")
+    c = pop_trips.evaluate("cost_per_person")
+    g = pop_trips.evaluate("ghg_per_person")
+    
+    v = pop_trips.evaluate("metrics_by_variable", variable="mode", plot=True)
+    v2 = pop_trips.evaluate("metrics_by_variable", variable="motive", plot=True)
+    v3 = pop_trips.evaluate("metrics_by_variable", variable="time_bin", plot=True)
+    v4 = pop_trips.evaluate("metrics_by_variable", variable="distance_bin", plot=True)
+    
diff --git a/mobility/choice_models/population_trips_parameters.py b/mobility/choice_models/population_trips_parameters.py
index 85e8425d..d6467310 100644
--- a/mobility/choice_models/population_trips_parameters.py
+++ b/mobility/choice_models/population_trips_parameters.py
@@ -1,24 +1,128 @@
-from dataclasses import dataclass
+from dataclasses import dataclass, fields
+from typing import List, Union
 
-@dataclass(frozen=True)
-class PopulationTripsParameters:
-    n_iterations: int = 1
+from mobility import Parameter, ParameterSet
+
+
+Number = Union[int, float]    
+
+
+@dataclass
+class PopulationTripsParameters(ParameterSet):
+    """To add a new parameter:
+        - add name, type and default value here,
+        - add the description in the __post_init__ method
+        - if the parameter is linked to another one, maybe add a constraint
+          in _validate_param_interdependency()
+        """
+    n_iterations: int = 4
     alpha: float = 0.01
     k_mode_sequences: int = 6
     dest_prob_cutoff: float = 0.99
     n_iter_per_cost_update: int = 3
     cost_uncertainty_sd: float = 1.0
     seed: int = 0
+    seeds: List[int] = None
     mode_sequence_search_parallel: bool = True
     min_activity_time_constant: float = 1.0
     simulate_weekend: bool = False
+    name: str = ""
+
+    
+    def __post_init__(self):
+        self.parameters = {
+    "param_n_iterations":         Parameter(name="Number of iterations", name_fr="Nombre d'itérations",
+                                     value=self.n_iterations,
+                                     description=("""Number of simulation iterations used to compute the population trips.
+                                                     Increase this to get more diverse programmes and to allow congestion feedbacks to propagate"""),
+                                     parameter_type=int, min_value=1, max_value=42, interval=1,
+                                     source_default="experience", parameter_role="simulation_parameter"),
+                                     
+    "param_alpha" :               Parameter(name="Alpha", name_fr="Alpha", description="Alpha",
+                                          value=self.alpha,
+                                          parameter_type=float, min_value=0.0, max_value=1.0, interval=0.005,
+                                          source_default="experience", parameter_role="simulation_parameter"),
+    "param_k_mode_sequences" :    Parameter(name="Number k of sampled mode sequences",
+                                          value=self.k_mode_sequences,
+                                          name_fr="Nombre k de séquences échantillonnées",
+                                          description="Increase this to get more mode sequence options for each destination sequence, decrease to speed up the simulation",
+                                          parameter_type=int, min_value=1, interval=1,
+                                          source_default="experience", parameter_role="simulation_parameter"),
+    "param_dest_prob_cutoff" :    Parameter(name="TBD", name_fr="TBD",
+                                          value=self.dest_prob_cutoff,
+                                          description="Increase this to get more destinations possibilities, decrease to speed up the simulation",
+                                          parameter_type=float, min_value=0.10, max_value=1.00, interval=0.05,
+                                          source_default="experience", parameter_role="simulation_parameter"),
+    "param_n_iter_per_cost_update" : Parameter(name="Number n of iteration between two cost updates",
+                                             value=self.n_iter_per_cost_update,
+                                             name_fr="Nombre n d'itérations entre deux mises à jours des coûts",
+                                             description="""Every n iterations, the model will update the costs:
+                                                         set to zero to ignore congestion in the simulation,
+                                                         set to 1 to update congestion at each iteration,
+                                                         set to a higher level to speed up the simulation""",
+                                             parameter_type=int, min_value=0, max_value=99, interval=1,
+                                             source_default="experience", parameter_role="simulation_parameter"),
+    "param_cost_uncertainty_sd" : Parameter(name="TBD", name_fr="TBD",
+                                          value=self.cost_uncertainty_sd,
+                                          description="Increase this to make the location of opportunities more uncertain around the destination points",
+                                          parameter_type=float, min_value=0.1, interval=0.1,
+                                          source_default="experience", parameter_role="simulation_parameter"),
+    "param_seed" :                Parameter(name="Seed", name_fr="Graine (seed)",
+                                          value=self.seed,
+                                          description="""Unique seed used for these trips.
+                                                      Change this value to get new programmes and results for a given set of inputs,
+                                                      keep this to a set value to get reproductible results.""",
+                                          parameter_type=int, min_value=0,
+                                          parameter_role="randomisation"),
+    "param_seeds" :               Parameter(name="Set of seeds", name_fr="Ensemble de graines (seeds)",
+                                          value=self.seeds,
+                                          description="Set of seeds, used to improve resultats by minimising noise",
+                                          parameter_type=list,
+                                          parameter_role="randomisation"),
+    "param_mode_sequence_search_parallel" : Parameter(name="TBD", name_fr="TBD",
+                                                    value=self.mode_sequence_search_parallel,
+                                                    description="Set to False to debug or for small simulations, otherwise set to True to speed up the simulation",
+                                                    parameter_type=bool,
+                                                    source_default="tests", parameter_role="simulation_time_parameter"),
+    "param_min_activity_time_constant" : Parameter(name="TBD", name_fr="TBD", description="TBD",
+                                                 value=self.min_activity_time_constant,
+                                                 parameter_type=float, interval=0.1,
+                                                 source_default="experience", parameter_role="simulation_parameter"),
+    "param_simulate_weekend" :    Parameter(name="Simulate weekend", name_fr="Simulation du week-end",
+                                          value=self.simulate_weekend,
+                                          description="Simulate or not trips for the week-ends",
+                                          parameter_type=bool,
+                                          source_default="speed_optimisation", parameter_role="simulation_parameter"),
+    "param_name": Parameter(name="Name", name_fr="Nom", description="Name of the run with those parameters",
+                            value=self.name, parameter_type=str)}
+ 
+    def to_dict(self):
+        # Return a structured description of all parameters
+        return {
+            "n_iterations": self.n_iterations,
+            "alpha": self.alpha,
+            "k_mode_sequences": self.k_mode_sequences,
+            "dest_prob_cutoff": self.dest_prob_cutoff,
+            "n_iter_per_cost_update": self.n_iter_per_cost_update,
+            "cost_uncertainty_sd": self.cost_uncertainty_sd,
+            "seed": self.seed,
+            "seeds": self.seeds,
+            "mode_sequence_search_parallel": self.mode_sequence_search_parallel,
+            "min_activity_time_constant": self.min_activity_time_constant,
+            "simulate_weekend": self.simulate_weekend,
+            "name": self.name,
+            # Add parameters in structured form
+            "parameters": {key: param.to_dict() for key, param in self.parameters.items()}
+        }                                        
+                                            
+                                            
+    def _validate_param_interdependency(self):
+        if self.seeds is not None:
+            for seed in self.seeds:
+                if (not isinstance(seed, int) or seed < 0):
+                    raise TypeError(f"Seed must be a positive int and not {seed}")
+        if self.n_iter_per_cost_update > self.n_iterations:
+            raise ValueError(f"n_iter_per_cost_update ({self.n_iter_per_cost_update}) should not be higher than the number of iterations ({self.n_iterations})")
+         
 
-    def validate(self) -> None:
-        assert self.n_iterations >= 1
-        assert 0.0 < self.dest_prob_cutoff <= 1.0
-        assert self.alpha >= 0.0
-        assert self.k_mode_sequences >= 1
-        assert self.n_iter_per_cost_update >= 0
-        assert self.cost_uncertainty_sd > 0.0
-        assert self.seed >= 0
-        assert self.min_activity_time_constant >= 0
+    
diff --git a/mobility/choice_models/results.py b/mobility/choice_models/results.py
index 5b968dd3..2fb3e829 100644
--- a/mobility/choice_models/results.py
+++ b/mobility/choice_models/results.py
@@ -3,6 +3,7 @@
 import polars as pl
 import numpy as np
 import plotly.express as px
+import plotly.graph_objects as go
 
 from typing import Literal
 from mobility.choice_models.evaluation.travel_costs_evaluation import TravelCostsEvaluation
@@ -10,6 +11,8 @@
 from mobility.choice_models.evaluation.routing_evaluation import RoutingEvaluation
 from mobility.choice_models.evaluation.public_transport_network_evaluation import PublicTransportNetworkEvaluation
 
+from mobility.parsers.work_home_flows import WorkHomeFlows_fr
+
 class Results:
     
     def __init__(
@@ -53,13 +56,15 @@ def __init__(
             "trip_count_by_demand_group": self.trip_count_by_demand_group,
             "distance_per_person": self.distance_per_person,
             "ghg_per_person": self.ghg_per_person,
+            #"ghg_emissions_per_trip_per_person": self.ghg_emissions_per_trip, #Test
             "time_per_person": self.time_per_person,
             "cost_per_person": self.cost_per_person,
             "immobility": self.immobility,
             "car_traffic": self.car_traffic,
             "travel_costs": self.travel_costs,
             "routing": self.routing,
-            "public_transport_network": self.public_transport_network
+            "public_transport_network": self.public_transport_network,
+            "ssi": self.ssi
         }
         
         
@@ -583,11 +588,13 @@ def metric_per_person(
             plot: bool = False,
             mask_outliers: bool = False,
             compare_with = None,
-            plot_delta = False
+            plot_delta = False,
+            labels = None
         ):
-        """        
-        Aggregate total value and value per person by demand group for this metric.
-        Metric can be : cost, time, distance, ghg
+        """
+        TODO
+        
+        Aggregate total travel time and time per person by demand group.
         
         Parameters
         ----------
@@ -738,7 +745,7 @@ def metric_per_person(
                 if mask_outliers:
                     tz["delta"] = self.mask_outliers(tz["delta"])
                 
-                self.plot_map(tz, "delta", color_continuous_scale="RdBu_r", color_continuous_midpoint=0)                
+                self.plot_map(tz, "delta", color_continuous_scale="RdBu_r", color_continuous_midpoint=0, labels=labels)                
 
         
         return metric_per_groups_and_transport_zones
@@ -768,6 +775,10 @@ def distance_per_person(self, *args, **kwargs):
     def ghg_per_person(self, *args, **kwargs):
         return self.metric_per_person("ghg_emissions_per_trip", *args, **kwargs)
     
+    def ghg_emissions_per_trip(self, *args, **kwargs):
+        """Test for compat"""
+        return self.metric_per_person("ghg_emissions_per_trip", *args, **kwargs)
+    
     def time_per_person(self, *args, **kwargs):
         """
         Aggregate total travel time and time per person by demand group.
@@ -792,8 +803,17 @@ def time_per_person(self, *args, **kwargs):
     
  
     def cost_per_person(self, *args, **kwargs):
-        """
-        Aggregate total travel cost and cost per person by demand group.
+        """    self,
+            weekday: bool = True,
+            plot: bool = False,
+            mask_outliers: bool = False,
+            compare_with = None,
+            plot_delta = False
+        ):
+        
+        TODO
+        
+        Aggregate total travel time and time per person by demand group.
         
         Parameters
         ----------
@@ -806,16 +826,17 @@ def cost_per_person(self, *args, **kwargs):
         -------
         pl.DataFrame
             Grouped by ['home_zone_id', 'csp', 'n_cars'] with:
-            - cost: sum(cost * n_persons)
+            - cost: sum(cost * n_persons) * 60.0 (minutes)
             - n_persons: group size
-            - time_per_person: cost / n_persons
+            - time_per_person: time / n_persons
         """
         return self.metric_per_person("cost", *args, **kwargs)
 
     
  
     def plot_map(self, tz, value: str = None, motive: str = None, plot_method: str = "browser",
-                 color_continuous_scale="Viridis", color_continuous_midpoint=None):
+                 color_continuous_scale="Viridis", color_continuous_midpoint=None,
+                 labels=None):
         """
         Render a Plotly choropleth for a transport-zone metric.
         
@@ -835,6 +856,8 @@ def plot_map(self, tz, value: str = None, motive: str = None, plot_method: str =
         logging.getLogger("kaleido").setLevel(logging.WARNING)
         #plot_method="png"
         
+        
+        
         fig = px.choropleth(
             tz.drop(columns="geometry"),
             geojson=json.loads(tz.to_json()),
@@ -848,6 +871,19 @@ def plot_map(self, tz, value: str = None, motive: str = None, plot_method: str =
             title=motive,
             subtitle=motive
         )
+        if labels is not None:
+            print(labels[["geometry"]])
+            labels = labels.to_crs(4326)
+            print(labels[["geometry"]])
+            for index, row in labels.iterrows():
+                if row["prominence"] == 1:
+                    fig.add_trace(go.Scattergeo(lon=[row["geometry"].centroid.x], lat=[row["geometry"].centroid.y], text=[row["local_admin_unit_name"]], mode="text", textposition="top center", textfont=dict(size=24, color="black"),))
+                elif row["prominence"] <3:
+                    fig.add_trace(go.Scattergeo(lon=[row["geometry"].centroid.x], lat=[row["geometry"].centroid.y], text=[row["local_admin_unit_name"]], mode="text", textposition="top center", textfont=dict(size=18, color="black"),))
+                else:
+                    fig.add_trace(go.Scattergeo(lon=[row["geometry"].centroid.x], lat=[row["geometry"].centroid.y], text=[row["local_admin_unit_name"]], mode="text", textposition="top center", textfont=dict(size=12, color="black"),))
+                    #plt.annotate(row["local_admin_unit_name"], (row["x"], row["y"]), size=size[0], ha="center", va="center", color=color)
+
         fig.update_geos(fitbounds="geojson", visible=False)
         fig.update_layout(margin=dict(l=0,r=0,t=0,b=0))
         fig.show(plot_method)
@@ -888,6 +924,82 @@ def routing(self, *args, **kwargs):
     def public_transport_network(self, *args, **kwargs):
         return PublicTransportNetworkEvaluation(self).get(*args, **kwargs)
          
+    def ssi(self, threshold=20, *args, **kwargs):
+        transport_zones_df = pl.from_pandas(self.transport_zones.get()[["local_admin_unit_id", "transport_zone_id"]])
+        od_flows = self.weekday_states_steps.collect().join(
+            transport_zones_df, left_on="from", right_on="transport_zone_id")
+        od_flows = od_flows.join(
+            transport_zones_df, left_on="to", right_on="transport_zone_id", suffix=('_to'))
+
+        od_flows_work = (
+            od_flows
+            .filter(pl.col("motive") == "work")
+            .group_by("local_admin_unit_id", "local_admin_unit_id_to")
+            .agg(pl.col("n_persons").sum())
+            .rename({"local_admin_unit_id" : "local_admin_unit_id_from"})
+            )
+
+        whf = pl.from_pandas(WorkHomeFlows_fr().get())
+        od_flows_work = od_flows_work.join(whf, on=["local_admin_unit_id_from","local_admin_unit_id_to"])
+
+        od_flows_work_above_threshold = (
+            od_flows_work.filter((pl.col("insee_flows") > threshold) | (pl.col("n_persons") > threshold))
+            )
+
+        # Classic SSI formula
+        ssi = (
+            od_flows_work_above_threshold
+            .select(
+                (
+                    2 * pl.min_horizontal("n_persons", "insee_flows").sum()
+                    / (pl.col("n_persons").sum() + pl.col("insee_flows").sum())
+                ).alias(f"ssi-{threshold}")
+            )
+        )
+
+
+        # Modified SSI formula from Liu & Yan 2020
+        # https://doi.org/10.1038/s41598-020-61613-y
+        # modified with threshold
+        # removes i -> i
+        ssi_liu_yan_2020 = (
+            od_flows_work_above_threshold
+            # exclude i = j and flows under threshold
+            .filter(pl.col("local_admin_unit_id_from") != pl.col("local_admin_unit_id_to"))
+            # compute local Sørensen term
+            .with_columns(
+                (
+                    2 * pl.min_horizontal("n_persons", "insee_flows")
+                    / (pl.col("n_persons") + pl.col("insee_flows"))
+                ).alias("local_ssi")
+            )
+            # mean
+            .select(
+                (pl.col("local_ssi").mean()).alias("ssi_liu_yan_2020")
+            )
+        )
+
+        # Other SSI from Yan et al. 2014
+        # https://doi.org/10.1098/rsif.2014.0834
+        # modified with threshold
+        ssi_yan_2014 = (
+            od_flows_work_above_threshold
+            # compute local Sørensen term
+            .with_columns(
+                (
+                    2 * pl.min_horizontal("n_persons", "insee_flows")
+                    / (pl.col("n_persons") + pl.col("insee_flows"))
+                ).alias("local_ssi")
+            )
+            # mean
+            .select(
+                (pl.col("local_ssi").mean()).alias("ssi_yan_2014")
+            )
+            )
+
+        ssis = pl.concat([ssi, ssi_liu_yan_2020, ssi_yan_2014], how="horizontal")
         
-        
+        return ssis
+
+
         
diff --git a/mobility/choice_models/state_updater.py b/mobility/choice_models/state_updater.py
index 46f0bff3..b7e7faeb 100644
--- a/mobility/choice_models/state_updater.py
+++ b/mobility/choice_models/state_updater.py
@@ -62,7 +62,7 @@ def get_new_states(
             motive_dur,
             iteration,
             motives,
-            parameters.min_activity_time_constant,
+            parameters["min_activity_time_constant"],
             tmp_folders
         )
         
@@ -73,7 +73,7 @@ def get_new_states(
             home_night_dur,
             home_motive.value_of_time_stay_home,
             stay_home_state,
-            parameters.min_activity_time_constant
+            parameters["min_activity_time_constant"]
         )
         
         transition_prob = self.get_transition_probabilities(current_states, possible_states_utility)
diff --git a/mobility/choice_models/top_k_mode_sequence_search.py b/mobility/choice_models/top_k_mode_sequence_search.py
index 2fb76a46..128b5120 100644
--- a/mobility/choice_models/top_k_mode_sequence_search.py
+++ b/mobility/choice_models/top_k_mode_sequence_search.py
@@ -107,12 +107,12 @@ def run(self, iteration, costs_aggregator, tmp_folders, parameters):
                 leg_modes[(from_, to_)].append(mode)
         
         
-        if parameters.mode_sequence_search_parallel is False:
+        if parameters["mode_sequence_search_parallel"] is False:
             
             logging.info("Finding probable mode sequences for the spatialized trip chains...")
             
             compute_subtour_mode_probabilities_serial(
-                parameters.k_mode_sequences,
+                parameters["k_mode_sequences"],
                 unique_location_chains,
                 costs,
                 leg_modes,
@@ -146,7 +146,7 @@ def run(self, iteration, costs_aggregator, tmp_folders, parameters):
                         "python",
                         "-u",
                         str(resources.files('mobility') / "transport_modes" / "compute_subtour_mode_probabilities.py"),
-                        "--k_sequences", str(parameters.k_mode_sequences),
+                        "--k_sequences", str(parameters["k_mode_sequences"]),
                         "--location_chains_path", str(location_chains_path),
                         "--costs_path", str(costs_path),
                         "--leg_modes_path", str(leg_modes_path),
diff --git a/mobility/experiments/gtfs_experiments.py b/mobility/experiments/gtfs_experiments.py
new file mode 100644
index 00000000..84cf3d00
--- /dev/null
+++ b/mobility/experiments/gtfs_experiments.py
@@ -0,0 +1,15 @@
+import gtfs_kit
+
+k = gtfs_kit.read_feed("https://app.mecatran.com/utw/ws/gtfsfeed/static/lio?apiKey=2b160d626f783808095373766f18714901325e45&type=gtfs_lio", "m")
+
+trips = k.get_trips()
+dates = k.get_dates()
+routes = k.get_routes()
+shapes = k.get_shapes(as_gdf=True)
+shapes.plot()
+
+t = gtfs_kit.routes.build_route_timetable(k, '452', ["20251020"])
+
+ts = k.compute_trip_stats(["452"])
+rs = gtfs_kit.routes.compute_route_stats(feed=k, dates=["20251020"], trip_stats_subset=ts)
+
diff --git a/mobility/experiments/multiple_seeds.py b/mobility/experiments/multiple_seeds.py
new file mode 100644
index 00000000..e89c9ff0
--- /dev/null
+++ b/mobility/experiments/multiple_seeds.py
@@ -0,0 +1,96 @@
+import os
+import dotenv
+import mobility
+import polars as pl
+
+dotenv.load_dotenv()
+
+
+mobility.set_params(
+    # package_data_folder_path=os.environ["MOBILITY_PACKAGE_DATA_FOLDER"],
+    # project_data_folder_path=os.environ["MOBILITY_PROJECT_DATA_FOLDER"]
+        package_data_folder_path="D:/mobility-data",
+        project_data_folder_path="D:/test-09",
+        debug=True
+)
+
+transport_zones = mobility.TransportZones("fr-09122", radius = 40, level_of_detail=1)
+tz = transport_zones.get()
+emp = mobility.EMPMobilitySurvey()
+pop = mobility.Population(transport_zones, sample_size = 10000)
+modes = [mobility.CarMode(transport_zones), mobility.WalkMode(transport_zones), mobility.BicycleMode(transport_zones)]
+surveys = [emp]
+motives = [mobility.HomeMotive(), mobility.WorkMotive(), mobility.OtherMotive(population=pop)]
+
+list_seeds = [0, 42, 64, 1000, 8888, 9999, 67654, 3333333333, 56563563208387, 7345678900000743]
+print(len(list_seeds))
+list_pop = []
+
+pop_trips_base = mobility.PopulationTrips(
+    pop,
+    modes,
+    motives,
+    surveys,
+    parameters=mobility.PopulationTripsParameters(n_iterations=6, seed=0, k_mode_sequences=6)
+    )
+labels=pop_trips_base.get_prominent_cities()
+
+
+for seed in list_seeds:
+    print("\n SEED:", seed, "\n")
+    pop_trips = mobility.PopulationTrips(
+        pop,
+        modes,
+        motives,
+        surveys,
+        parameters=mobility.PopulationTripsParameters(n_iterations=6, seed=seed, k_mode_sequences=6)
+        )
+    cost_per_person = pop_trips.evaluate("cost_per_person", plot_delta=False, compare_with=pop_trips_base, labels=labels)
+    #dist_per_person = pop_trips.evaluate("distance_per_person", plot_delta=True, compare_with=pop_trips_base)
+    #time_per_person = pop_trips.evaluate("time_per_person", plot_delta=True, compare_with=pop_trips_base)
+    #ghg_per_person = pop_trips.evaluate("ghg_per_person", plot_delta=True, compare_with=pop_trips_base)
+    list_pop.append(pop_trips)
+    
+    
+cost_all = list_pop[0].evaluate("cost_per_person").select("transport_zone_id", "cost", "n_persons").group_by("transport_zone_id").agg(pl.col("cost").sum(), pl.col("n_persons").sum())
+i = 1
+for pop in list_pop[1:10]:
+    cost_add = pop.evaluate("cost_per_person").select("transport_zone_id", "cost", "n_persons").group_by("transport_zone_id").agg(pl.col("cost").sum(), pl.col("n_persons").sum())
+    suffix = "_" + str(i)
+    cost_all = cost_all.join(cost_add, on="transport_zone_id", suffix=suffix)
+    i += 1
+    
+cost_all = cost_all.sort(by="cost")
+
+pandas_cost=cost_all.select(pl.col("^cost.*$")).to_pandas()
+pandas_cost.plot()
+
+
+distance_all = list_pop[0].evaluate("distance_per_person").select("transport_zone_id", "distance", "n_persons").group_by("transport_zone_id").agg(pl.col("distance").sum(), pl.col("n_persons").sum())
+i = 1
+for pop in list_pop[1:10]:
+    distance_add = pop.evaluate("distance_per_person").select("transport_zone_id", "distance", "n_persons").group_by("transport_zone_id").agg(pl.col("distance").sum(), pl.col("n_persons").sum())
+    suffix = "_" + str(i)
+    distance_all = distance_all.join(distance_add, on="transport_zone_id", suffix=suffix)
+    i += 1
+    
+distance_all = distance_all.sort(by="distance")
+
+pandas_distance=distance_all.select(pl.col("^distance.*$")).to_pandas()
+pandas_distance.plot()
+
+time_all = list_pop[0].evaluate("time_per_person").select("transport_zone_id", "time", "n_persons").group_by("transport_zone_id").agg(pl.col("time").sum(), pl.col("n_persons").sum())
+i = 1
+for pop in list_pop[1:10]:
+    time_add = pop.evaluate("time_per_person").select("transport_zone_id", "time", "n_persons").group_by("transport_zone_id").agg(pl.col("time").sum(), pl.col("n_persons").sum())
+    suffix = "_" + str(i)
+    time_all = time_all.join(time_add, on="transport_zone_id", suffix=suffix)
+    i += 1
+    
+time_all = time_all.sort(by="time")
+
+pandas_time=time_all.select(pl.col("^time.*$")).to_pandas()
+pandas_time.plot()
+
+for pop in list_pop:
+    pop.plot_od_flows()
\ No newline at end of file
diff --git a/mobility/experiments/sensitivity.py b/mobility/experiments/sensitivity.py
new file mode 100644
index 00000000..b0f0d316
--- /dev/null
+++ b/mobility/experiments/sensitivity.py
@@ -0,0 +1,78 @@
+import os
+import dotenv
+import logging
+
+import mobility
+
+import polars as pl
+
+dotenv.load_dotenv()
+
+os.environ["MOBILITY_GTFS_DOWNLOAD_DATE"] = "2025/12/01"
+
+mobility.set_params(
+        package_data_folder_path="D:/mobility-data",
+        project_data_folder_path="D:/test-09",
+        debug=False
+)
+
+global_metrics = pl.DataFrame()
+ssis = pl.DataFrame()
+ssis200 = pl.DataFrame()
+
+for radius in range(36, 51, 2):
+    if radius not in [36]: #blockage for this one
+        print(f"\nRADIUS: {radius}")
+        transport_zones = mobility.TransportZones("fr-87085", radius = radius, level_of_detail=1)
+        emp = mobility.EMPMobilitySurvey()
+        pop = mobility.Population(transport_zones, sample_size = 1000)
+        modes = [mobility.CarMode(transport_zones), mobility.WalkMode(transport_zones),
+                 mobility.BicycleMode(transport_zones), mobility.PublicTransportMode(transport_zones)]
+        surveys = [emp]
+        motives = [mobility.HomeMotive(), mobility.WorkMotive(), mobility.OtherMotive(population=pop)]
+        
+        
+        # Simulating the trips for this population for three modes : car, walk and bicyle, and only home and work motives (OtherMotive is mandatory)
+        pop_trips = mobility.PopulationTrips(
+            pop,
+            modes,
+            motives,
+            surveys,
+            parameters=mobility.PopulationTripsParameters(n_iterations=4, k_mode_sequences=42)
+            )
+    
+        labels=pop_trips.get_prominent_cities()
+        
+        if global_metrics.is_empty():
+            global_metrics = pop_trips.evaluate("global_metrics")
+        else:
+            suffix = str(radius)
+            global_metrics = global_metrics.join(pop_trips.evaluate("global_metrics"), on =["country", "variable"], suffix=suffix)
+        
+        if radius % 40 == 0:            
+            metrics_by_mode = pop_trips.evaluate("metrics_by_variable", variable="mode", plot=True)
+            metrics_by_motive = pop_trips.evaluate("metrics_by_variable", variable="motive", plot=True)
+        
+        # # OD flows between transport zones
+        pop_trips.plot_od_flows(mode="car", level_of_detail=1, labels=labels)
+        pop_trips.plot_od_flows(mode="walk", level_of_detail=1, labels=labels)
+        pop_trips.plot_od_flows(mode="bicycle", level_of_detail=1, labels=labels)
+        pop_trips.plot_od_flows(mode="public_transport", labels=labels)
+        
+        if ssis.is_empty():
+            rad = pl.DataFrame({"radius": radius})
+            ssis = pl.concat([rad, pop_trips.evaluate("ssi")], how="horizontal")
+            ssis200 = pl.concat([rad, pop_trips.evaluate("ssi", threshold=200)], how="horizontal")
+        else:
+            rad = pl.DataFrame({"radius": radius})
+            ssis = pl.concat([ssis, pl.concat([rad, pop_trips.evaluate("ssi")], how="horizontal")])
+            ssis200 = pl.concat([ssis200, pl.concat([rad, pop_trips.evaluate("ssi", threshold=200)], how="horizontal")])
+        
+        print(global_metrics)
+        print(ssis)
+        print(ssis200)
+ssis.to_pandas().plot(x="radius")
+ssis200.to_pandas().plot(x="radius")
+    
+    
+        
\ No newline at end of file
diff --git a/mobility/choice_models/destination_choice_model.py b/mobility/legacy/destination_choice_model.py
similarity index 100%
rename from mobility/choice_models/destination_choice_model.py
rename to mobility/legacy/destination_choice_model.py
diff --git a/mobility/choice_models/leisure_destination_choice_model.py b/mobility/legacy/leisure_destination_choice_model.py
similarity index 100%
rename from mobility/choice_models/leisure_destination_choice_model.py
rename to mobility/legacy/leisure_destination_choice_model.py
diff --git a/mobility/choice_models/shopping_destination_choice_model.py b/mobility/legacy/shopping_destination_choice_model.py
similarity index 100%
rename from mobility/choice_models/shopping_destination_choice_model.py
rename to mobility/legacy/shopping_destination_choice_model.py
diff --git a/mobility/choice_models/studies_destination_choice_model.py b/mobility/legacy/studies_destination_choice_model.py
similarity index 100%
rename from mobility/choice_models/studies_destination_choice_model.py
rename to mobility/legacy/studies_destination_choice_model.py
diff --git a/mobility/model_parameters.py b/mobility/model_parameters.py
new file mode 100644
index 00000000..3afba3e5
--- /dev/null
+++ b/mobility/model_parameters.py
@@ -0,0 +1,115 @@
+from dataclasses import dataclass, field, fields, asdict
+from typing import List, Union
+
+Number = Union[int, float]    
+
+@dataclass(frozen=True)
+class Parameter:
+    name: str
+    name_fr: str
+    value: Number | bool
+    description: str
+    parameter_type: type
+    #default_value: Number | bool
+    possible_values: List[float] | List[int] | tuple = None
+    min_value: Number = None
+    max_value: Number = None
+    unit: str = None
+    interval: Number = None
+    source_default: str = ""
+    parameter_role: str = ""
+    
+    def to_dict(self):
+        # Convert the parameter to a dictionary with serializable values
+        return {
+            "name": self.name,
+            "name_fr": self.name_fr,
+            "value": self.value,
+            "description": self.description,
+            "parameter_type": str(self.parameter_type),  # Convert type to string
+            "possible_values": self.possible_values,
+            "min_value": self.min_value,
+            "max_value": self.max_value,
+            "unit": self.unit,
+            "interval": self.interval,
+            "source_default": self.source_default,
+            "parameter_role": self.parameter_role,
+        }
+
+    # def get(self):
+    #     """Return parameter value."""
+    #     val = self.default_value
+    #     self._validate(val)
+    #     return val
+
+    def set(self, new_value):
+        self.value = new_value
+
+    def validate(self):
+        if self.value is None: #todo: improve this!
+            return None
+        if self.parameter_type is not None:
+            if not isinstance(self.value, self.parameter_type):
+                t = type(self.value)
+                raise TypeError(f"Parameter '{self.name}' must be {self.parameter_type} (currently {t}).")
+
+        if self.min_value is not None and self.value < self.min_value:
+            raise ValueError(
+                f"Parameter '{self.name}' below minimum {self.min_value}"
+            )
+
+        if self.max_value is not None and self.value > self.max_value:
+            raise ValueError(
+                f"Parameter '{self.name}' above maximum {self.max_value}"
+            )
+
+    def __repr__(self):
+        unit_str = f" [{self.unit}]" if self.unit else ""
+        return f"<Parameter {self.name}={self.value}{unit_str}>"
+
+
+    def get_values_for_sensitivity_analyses(self, i_max=10):
+        value = self.value
+        values = [value]
+        if self.interval is None:
+            raise ValueError("To run a sensitivity analysis, interval must be specified in the parameter configuration")
+        i = 0
+        if self.max_value is not None:
+            while i < i_max and value < self.max_value:
+                value += self.interval
+                values.append(round(value,3))
+                i += 1
+        else:
+            while i < i_max:
+                value += self.interval
+                values.append(round(value,3))                
+                i += 1
+        value = self.value
+        i = 0
+        
+        if self.min_value is not None:
+            while i < i_max and value > self.min_value:
+                value -= self.interval
+                values.append(round(value,3))
+                i += 1
+        else:
+            while i < i_max:
+                value -= self.interval
+                values.append(round(value,3))                
+                i += 1
+                
+        values.sort()
+        return values
+    
+@dataclass    
+class ParameterSet:
+    parameters : dict = field(init=False, compare=False)    
+    
+    def validate(self):
+        for param in fields(self)[1:]:
+            param_name = "param_" + param.name
+            self.parameters[param_name].validate()
+        self._validate_param_interdependency()
+            
+    def _validate_param_interdependency(self):
+        pass
diff --git a/mobility/motives/other.py b/mobility/motives/other.py
index 7e461d15..cc05c349 100644
--- a/mobility/motives/other.py
+++ b/mobility/motives/other.py
@@ -11,7 +11,7 @@ def __init__(
         value_of_time: float = 10.0,
         saturation_fun_ref_level: float = 1.5,
         saturation_fun_beta: float = 4.0,
-        radiation_lambda: float = 0.99986,
+        radiation_lambda: float = 0.9999999,
         opportunities: pd.DataFrame = None,
         population: Population = None
     ):
diff --git a/mobility/motives/work.py b/mobility/motives/work.py
index e8834c8d..4b400b77 100644
--- a/mobility/motives/work.py
+++ b/mobility/motives/work.py
@@ -14,7 +14,7 @@ def __init__(
         saturation_fun_ref_level: float = 1.5,
         saturation_fun_beta: float = 4.0,
         survey_ids: List[str] = ["9.91"],
-        radiation_lambda: float = 0.99986,
+        radiation_lambda: float = 0.9999999,
         opportunities: pd.DataFrame = None,
         utilities: pd.DataFrame = None,
         country_utilities: Dict = None
diff --git a/mobility/parsers/__init__.py b/mobility/parsers/__init__.py
index d94ae8d7..7645625b 100644
--- a/mobility/parsers/__init__.py
+++ b/mobility/parsers/__init__.py
@@ -1,10 +1,11 @@
 from .patch_openpyxl import patch_openpyxl
 patch_openpyxl()
 
-from .work_home_flows import download_work_home_flows
+#from .work_home_flows import download_work_home_flows
 from .city_legal_population import CityLegalPopulation
 from .census_localized_individuals import CensusLocalizedIndividuals
 from .mobility_survey import MobilitySurvey
 from .local_admin_units import LocalAdminUnits
 from .jobs_active_population_distribution import JobsActivePopulationDistribution
-from .jobs_active_population_flows import JobsActivePopulationFlows
\ No newline at end of file
+from .jobs_active_population_flows import JobsActivePopulationFlows
+from .work_home_flows import WorkHomeFlows_fr
\ No newline at end of file
diff --git a/mobility/parsers/work_home_flows.py b/mobility/parsers/work_home_flows.py
index affc6c4d..4200329f 100644
--- a/mobility/parsers/work_home_flows.py
+++ b/mobility/parsers/work_home_flows.py
@@ -1,38 +1,106 @@
+import logging
+import numpy as np
 import pandas as pd
 import os
 from pathlib import Path
+import pathlib
+
 import requests
 import zipfile
 
+from mobility.file_asset import FileAsset
+from mobility.parsers.download_file import download_file
 
-def download_work_home_flows(proxies={}):
-    """
-    Downloads (if needed) the INSEE census mobility flows data
-    (hosted on data.gouv.fr by the mobility project at
-    https://www.data.gouv.fr/fr/datasets/r/f3f22487-22d0-45f4-b250-af36fc56ccd0)
-    """
 
-    data_folder_path = (
-        Path(os.path.dirname(__file__)).parents[0] / "data/insee/work_home_flows"
-    )
+#script
+import mobility
 
-    if data_folder_path.exists() is False:
-        os.makedirs(data_folder_path)
 
-    # Download the raw survey data from insee.fr if needed
-    path = data_folder_path / "rp2019-mobpro-csv.zip"
 
-    if path.exists() is False:
-        # Download the zip file
-        r = requests.get(
-            url="https://www.data.gouv.fr/fr/datasets/r/f3f22487-22d0-45f4-b250-af36fc56ccd0",
-            proxies=proxies,
-        )
-        with open(path, "wb") as file:
-            file.write(r.content)
+class WorkHomeFlows_fr(FileAsset):
+
+    
+    def __init__(self, year="2021"):
+        
+        inputs = {"year": year}
+
+        file_name = f"insee_mobpro_{year}.parquet"
+        cache_path = pathlib.Path(os.environ["MOBILITY_PACKAGE_DATA_FOLDER"]) / "insee" / "flows" / file_name
 
-        # Unzip the content
+        super().__init__(inputs, cache_path)
+        
+    def get_cached_asset(self) -> pd.DataFrame:
+
+        logging.info("French home-work flows already prepared. Reusing the file : " + str(self.cache_path))
+        flows = pd.read_parquet(self.cache_path)
+
+        return flows
+    
+    
+    def create_and_get_asset(self) -> pd.DataFrame:
+        """
+        Parse and format french home-work flows for the given year.
+        
+        Returns:
+            A pandas.DataFrame giving the french home-work flows
+        """
+
+        urls ={
+        "2021" : "https://www.insee.fr/fr/statistiques/fichier/8201899/base-flux-mobilite-domicile-lieu-travail-2021-csv.zip",
+        "2020" : "https://www.insee.fr/fr/statistiques/fichier/7630376/base-flux-mobilite-domicile-lieu-travail-2020-csv.zip",
+        "2019" : "https://www.insee.fr/fr/statistiques/fichier/6454112/base-csv-flux-mobilite-domicile-lieu-travail-2019.zip",
+        "2018" : "https://www.insee.fr/fr/statistiques/fichier/5393835/base-csv-flux-mobilite-domicile-lieu-travail-2018.zip"
+        }
+        
+        year = self.year
+        
+        folder = pathlib.Path(os.environ["MOBILITY_PACKAGE_DATA_FOLDER"]) / "insee" / "flows"
+        if folder.exists() is False:
+            os.mkdir(folder)
+        path = folder / f"insee_mobpro_{year}.zip"
+        download_file(urls[year], path)
+        
         with zipfile.ZipFile(path, "r") as zip_ref:
-            zip_ref.extractall(data_folder_path)
+            zip_ref.extractall(folder)
+            
+        match self.year:
+            case "2021":
+                file_name= "base-flux-mobilite-domicile-lieu-travail-2021.csv"
+                col_name = "NBFLUX_C21_ACTOCC15P"
+            case "2020":
+                file_name= "base-flux-mobilite-domicile-lieu-travail-2020.csv"
+                col_name = "NBFLUX_C20_ACTOCC15P"
+            case "2019":
+                file_name= "base-flux-mobilite-domicile-lieu-travail-2019.csv"
+                col_name = "NBFLUX_C19_ACTOCC15P"
+            case "2018":
+                file_name= "base-flux-mobilite-domicile-lieu-travail-2018.csv"
+                col_name = "NBFLUX_C18_ACTOCC15P"                
+                
+        flows = pd.read_csv(
+            folder / file_name,
+            sep=";",
+            usecols=["CODGEO", "DCLT", col_name],
+            dtype={"CODGEO": str, "DCLT": str, col_name: np.float32}
+        )
+        flows.columns = ["local_admin_unit_id_from", "local_admin_unit_id_to","insee_flows"]
+        
+        flows["local_admin_unit_id_from"] = "fr-" + flows["local_admin_unit_id_from"]
+        flows["local_admin_unit_id_to"] = "fr-" + flows["local_admin_unit_id_to"]
+ 
+        flows.to_parquet(self.cache_path)
+        
+        return flows
 
-    return str(data_folder_path / "FD_MOBPRO_2019.csv")
+if __name__ == "__main__":
+    
+    
+    mobility.set_params(
+        # package_data_folder_path=os.environ["MOBILITY_PACKAGE_DATA_FOLDER"],
+        # project_data_folder_path=os.environ["MOBILITY_PROJECT_DATA_FOLDER"]
+            package_data_folder_path="D:/mobility-data",
+            project_data_folder_path="D:/test-09",
+            debug=True
+    )
+    f = WorkHomeFlows_fr(year="2018")
+    f1 = f.get()
\ No newline at end of file
diff --git a/mobility/transport_modes/bicycle/bicycle_mode.py b/mobility/transport_modes/bicycle/bicycle_mode.py
index 888fed22..b8ba332f 100644
--- a/mobility/transport_modes/bicycle/bicycle_mode.py
+++ b/mobility/transport_modes/bicycle/bicycle_mode.py
@@ -35,8 +35,8 @@ def __init__(
         
         if generalized_cost_parameters is None:
             generalized_cost_parameters = GeneralizedCostParameters(
-                cost_constant=0.0,
-                cost_of_distance=0.0,
+                cost_constant=25.0,
+                cost_of_distance=0.20,
                 cost_of_time=CostOfTimeParameters()
             )
         
diff --git a/mobility/transport_modes/car/car_mode.py b/mobility/transport_modes/car/car_mode.py
index 9f51e8ed..3f868066 100644
--- a/mobility/transport_modes/car/car_mode.py
+++ b/mobility/transport_modes/car/car_mode.py
@@ -45,8 +45,8 @@ def __init__(
         
         if generalized_cost_parameters is None:
             generalized_cost_parameters = GeneralizedCostParameters(
-                cost_constant=0.0,
-                cost_of_distance=0.1,
+                cost_constant=2.0,
+                cost_of_distance=0.05,
                 cost_of_time=CostOfTimeParameters()
             )
             
diff --git a/mobility/transport_modes/modal_transfer.py b/mobility/transport_modes/modal_transfer.py
index 9568a95b..a13cfd58 100644
--- a/mobility/transport_modes/modal_transfer.py
+++ b/mobility/transport_modes/modal_transfer.py
@@ -16,11 +16,11 @@ class IntermodalTransfer():
     
     # Max travel time from or to the connection nodes, estimated from the crow
     # fly distance around connection nodes and an average speed
-    max_travel_time: float # hours
-    average_speed: float # km/h
+    max_travel_time: float = 20.0 / 60.0 # hours
+    average_speed: float = 4.0 # km/h
     
     # Average transfer time between the two connected modes
-    transfer_time: float # minutes
+    transfer_time: float = 1.0 # minutes
     
     # Optional shortcuts to make some connections faster than average
     shortcuts_transfer_time: float = None # minutes
diff --git a/mobility/transport_modes/public_transport/gtfs/gtfs_data.py b/mobility/transport_modes/public_transport/gtfs/gtfs_data.py
index 93f9f4c4..b2b85b12 100644
--- a/mobility/transport_modes/public_transport/gtfs/gtfs_data.py
+++ b/mobility/transport_modes/public_transport/gtfs/gtfs_data.py
@@ -63,13 +63,20 @@ def is_gtfs_file_ok(self, path):
         if "e8f2aceaaaa2493f6041dc7f0251f325-5d7ae44c16ad373ca1afbc4590f53256_gtfs-2015-chamonix-mobilit" in path.name:
             logging.info("Manual exception, GTFS not used from path", path)
             return False
+        if "datasud" in path.name:
+            logging.info("Manual exception, GTFS not used from path", path)
+            return False
+        if "4a590cb87669fe2bc39328652ef1d2e9_gtfs_generic_eu" in path.name:
+            logging.info(f"Manual exception, Flixbus GTFS not used from path {path.name}")
+            return False
+
         
         try:
             with zipfile.ZipFile(path, 'r') as zip_ref:
                 zip_contents = zip_ref.namelist()
             has_an_agency = "agency.txt" in zip_contents
             if has_an_agency:
-                logging.info("Downloaded file is a proper GTFS zip which contains an agency file.")
+                logging.debug("Downloaded file is a proper GTFS zip which contains an agency file.")
             else:
                 logging.info("Downloaded file is a proper GTFS zip but does not contain an agency file, it will not be used by Mobility.")
             return True
diff --git a/mobility/transport_modes/public_transport/public_transport_mode.py b/mobility/transport_modes/public_transport/public_transport_mode.py
index 2837f046..b695b770 100644
--- a/mobility/transport_modes/public_transport/public_transport_mode.py
+++ b/mobility/transport_modes/public_transport/public_transport_mode.py
@@ -4,6 +4,7 @@
 
 from mobility.transport_zones import TransportZones
 from mobility.transport_modes.transport_mode import TransportMode
+from mobility.transport_modes.walk import WalkMode
 from mobility.transport_modes.public_transport.public_transport_routing_parameters import PublicTransportRoutingParameters
 from mobility.transport_modes.public_transport.public_transport_travel_costs import PublicTransportTravelCosts
 from mobility.transport_modes.public_transport.public_transport_generalized_cost import PublicTransportGeneralizedCost
@@ -32,8 +33,8 @@ def __init__(
         transport_zones: TransportZones,
         first_leg_mode: TransportMode = None,
         last_leg_mode: TransportMode = None,
-        first_intermodal_transfer: IntermodalTransfer = None,
-        last_intermodal_transfer: IntermodalTransfer = None,
+        first_intermodal_transfer: IntermodalTransfer = IntermodalTransfer(),
+        last_intermodal_transfer: IntermodalTransfer = IntermodalTransfer(),
         routing_parameters: PublicTransportRoutingParameters = PublicTransportRoutingParameters(),
         generalized_cost_parameters: GeneralizedCostParameters = None,
         survey_ids: List[str] = [
@@ -44,6 +45,11 @@ def __init__(
         ghg_intensity: float = 0.05
     ):
         
+        if first_leg_mode is None:
+            first_leg_mode = WalkMode(transport_zones)
+        if last_leg_mode is None:
+            last_leg_mode = WalkMode(transport_zones)
+            
         travel_costs = PublicTransportTravelCosts(
             transport_zones,
             routing_parameters,
@@ -57,8 +63,8 @@ def __init__(
         
         if generalized_cost_parameters is None:
             generalized_cost_parameters = GeneralizedCostParameters(
-                cost_constant=0.0,
-                cost_of_distance=0.1,
+                cost_constant=1.0,
+                cost_of_distance=0.0,
                 cost_of_time=CostOfTimeParameters()
             )
         
diff --git a/mobility/transport_modes/walk/walk_mode.py b/mobility/transport_modes/walk/walk_mode.py
index 7a936603..a954c6f2 100644
--- a/mobility/transport_modes/walk/walk_mode.py
+++ b/mobility/transport_modes/walk/walk_mode.py
@@ -34,8 +34,8 @@ def __init__(
         
         if generalized_cost_parameters is None:
             generalized_cost_parameters = GeneralizedCostParameters(
-                cost_constant=0.0,
-                cost_of_distance=0.0,
+                cost_constant=0.5,
+                cost_of_distance=0.1,
                 cost_of_time=CostOfTimeParameters()
             )