diff --git a/docs/optimization.md b/docs/optimization.md
index c343d1e..9dae76a 100644
--- a/docs/optimization.md
+++ b/docs/optimization.md
@@ -4,22 +4,23 @@
 
 ### Log posterior probability
 
-***class* log_posterior_probability(model, parameter_names, bounds, data, states, log_likelihood_fnc, log_likelihood_fnc_args, weights=None, log_prior_prob_fnc=None, log_prior_prob_fnc_args=None, initial_states=None, aggregation_function=None, labels=None)**
+***class* log_posterior_probability(model, parameter_names, bounds, data, states, log_likelihood_fnc, log_likelihood_fnc_args, start_sim=None, weights=None, log_prior_prob_fnc=None, log_prior_prob_fnc_args=None, initial_states=None, aggregation_function=None, labels=None)**
 
 **Parameters:**
 
-* **model** (object) - An initialized ODE or JumpProcess.
+* **model** (object) - An initialized `pySODM.models.base.ODE` or `pySODM.models.base.JumpProcess` model.
 * **parameter_names** (list) - Names of model parameters (type: str) to calibrate. Model parameters must be of type float (0D), list containing float (1D), or np.ndarray (nD).
-* **bounds** (list) - Lower and upper bound of calibrated parameters provided as lists/tuples containing lower and upper bound: example: `[(lb_1, ub_1), ..., (lb_n, ub_n)]`. Values falling outside these bounds will be restricted to the provided ranges before simulating the model.
-* **data** (list) - Contains the datasets (type: pd.Series/pd.DataFrame) the model should be calibrated to. For one dataset use `[dataset,]`. Must contain a time index named `time` or `date`. Additional axes must be implemented using a `pd.Multiindex` and must bear the names/contain the coordinates of a valid model dimension.
-* **states** (list) - Names of the model states (type: str) the respective datasets should be matched to.
-* **log_likelihood_fnc** (list) - Contains a log likelihood function for every provided dataset.
-* **log_likelihood_fnc_args** (list) - Contains the arguments of the log likelihood functions. If the log likelihood function has no arguments (`ll_poisson`), provide an empty list.
+* **bounds** (list) - Lower and upper bound of calibrated parameters. Provided as a list or tuple containing lower and upper bound: example: `bounds = [(lb_1, ub_1), ..., (lb_n, ub_n)]`.
+* **data** (list) - Contains the datasets (type: pd.Series/pd.DataFrame) the model should be calibrated to. If there is only one dataset use `data = [df,]`. Dataframe must contain an index named `time` or `date`. Stratified data can be incorporated using a `pd.Multiindex`, whose index levels must have names corresponding to valid model dimensions, and whose indices must be valid dimension coordinates.
+* **states** (list) - Names of the model states (type: str) the respective datasets should be matched to. Must have the same length as `data`.
+* **log_likelihood_fnc** (list) - Contains a log likelihood function for every provided dataset. Must have the same length as `data`.
+* **log_likelihood_fnc_args** (list) - Contains the arguments of the log likelihood functions. If the log likelihood function has no arguments (such as `ll_poisson`), provide an empty list. Must have the same length as `data`.
+* **start_sim** (int/float or str/datetime) - optional - Can be used to alter the start of the simulation. By default, the start of the simulation is chosen as the earliest time/date found in the datasets. 
 * **weights** (list) - optional - Contains the weights of every dataset in the final log posterior probability. Defaults to one for every dataset.
 * **log_prior_prob_fnc** (list) - optional - Contains a prior probability function for every calibrated parameter. Defaults to a uniform prior using the provided bounds.
 * **log_prior_prob_fnc_args** (list) - optional - Contains the arguments of the provided prior probability functions.
-* **initial_states** (list) - optional - Contains a dictionary of initial states for every dataset. 
-* **aggregation_function** (callable function or list) - optional - A user-defined function to manipulate the model output before matching it to data. The function takes as input an `xarray.DataArray`, resulting from selecting the simulation output at the state we wish to match to the dataset (`model_output_xarray_Dataset['state_name']`), as its input. The output of the function must also be an `xarray.DataArray`. No checks are performed on the input or output of the aggregation function, use at your own risk. Illustrative use case: I have a spatially explicit epidemiological model and I desire to simulate it a high spatial resolutioni. However, data is only available on a lower level of spatial resolution. Hence, I use an aggregation function to properly aggregate the spatial levels. I change the coordinates on the spatial dimensions in the model output. Valid inputs for the argument `aggregation_function`are: 1) one callable function --> applied to every dataset. 2) A list containing one callable function --> applied to every dataset. 3) A list containing a callable function for every dataset --> every dataset has its own aggregation function.
+* **initial_states** (list) - optional - Contains a dictionary of initial states for every dataset.
+* **aggregation_function** (callable function or list) - optional - A user-defined function to manipulate the model output before matching it to data. The function takes as input an `xarray.DataArray`, resulting from selecting the simulation output at the state we wish to match to the dataset (`model_output_xarray_Dataset['state_name']`), as its input. The output of the function must also be an `xarray.DataArray`. No checks are performed on the input or output of the aggregation function, use at your own risk. Illustrative use case: I have a spatially explicit epidemiological model and I desire to simulate it a fine spatial resolution. However, data is only available on a coarser level. Hence, I use an aggregation function to properly aggregate the spatial levels. I change the coordinates on the spatial dimensions in the model output. Valid inputs for the argument `aggregation_function`are: 1) one callable function --> applied to every dataset. 2) A list containing one callable function --> applied to every dataset. 3) A list containing a callable function for every dataset --> every dataset has its own aggregation function.
 * **labels** (list) - optional - Contains a custom label for the calibrated parameters. Defaults to the names provided in `parameter_names`.
 
 **Methods:**
diff --git a/docs/quickstart.md b/docs/quickstart.md
index cec923c..6040ae6 100644
--- a/docs/quickstart.md
+++ b/docs/quickstart.md
@@ -47,15 +47,19 @@ To initialize the model, provide a dictionary containing the initial values of t
 model = SIR(states={'S': 1000, 'I': 1}, parameters={'beta': 0.35, 'gamma': 5})
 ```
 
-Simulate the model using its `sim()` method. pySODM supports the use of `datetime` types as timesteps.
+Simulate the model using its `sim()` method. pySODM supports the use of dates to index simulations, string representations of dates with the format `'yyyy-mm-dd'` as well as `datetime.datetime()` can be used. 
 
 ```python
 # Timesteps
 out = model.sim(121)
 
-# Dates
+# String representation of dates: 'yyyy-mm-dd' only
 out = model.sim(['2022-12-01', '2023-05-01'])
 
+# Datetime representation of time + date
+from datetime import datetime as datetime
+out = model.sim([datetime(2022, 12, 1), datetime(2023, 5, 1)])
+
 # Tailor method and tolerance of integrator
 out = model.sim(121, method='RK45', rtol='1e-4')
 ```
diff --git a/src/pySODM/models/base.py b/src/pySODM/models/base.py
index adf24e6..4b34507 100644
--- a/src/pySODM/models/base.py
+++ b/src/pySODM/models/base.py
@@ -11,7 +11,7 @@
 from pySODM.models.validation import merge_parameter_names_parameter_stratified_names, validate_draw_function, validate_simulation_time, validate_dimensions, \
                                         validate_time_dependent_parameters, validate_integrate, check_duplicates, build_state_sizes_dimensions, validate_dimensions_per_state, \
                                             validate_initial_states, validate_integrate_or_compute_rates_signature, validate_provided_parameters, validate_parameter_stratified_sizes, \
-                                                validate_apply_transitionings_signature, validate_compute_rates, validate_apply_transitionings
+                                                validate_apply_transitionings_signature, validate_compute_rates, validate_apply_transitionings, validate_solution_methods_ODE, validate_solution_methods_JumpProcess
 
 class JumpProcess:
     """
@@ -375,7 +375,7 @@ def _mp_sim_single(self, drawn_parameters, seed, time, actual_start_date, method
         np.random.seed(seed)
         return self._sim_single(time, actual_start_date, method, tau, output_timestep)
 
-    def sim(self, time, warmup=0, N=1, draw_function=None, draw_function_kwargs={}, processes=None, method='tau_leap', tau=0.5, output_timestep=1):
+    def sim(self, time, warmup=0, N=1, draw_function=None, draw_function_kwargs={}, processes=None, method='tau_leap', tau=1, output_timestep=1):
 
         """
         Simulate a model during a given time period.
@@ -422,27 +422,15 @@ def sim(self, time, warmup=0, N=1, draw_function=None, draw_function_kwargs={},
         output: xarray.Dataset
             Simulation output
         """
-
-        # Input checks on solution method and timestep
-        if not isinstance(method, str):
-            raise TypeError(
-                "solver method 'method' must be of type string"
-                )
-        if not isinstance(tau, (int,float)):
-            raise TypeError(
-                "discrete timestep 'tau' must be of type int or float"
-                )
+        
+        # Input checks on solution settings
+        validate_solution_methods_JumpProcess(method, tau)
         # Input checks on supplied simulation time
         time, actual_start_date = validate_simulation_time(time, warmup)
         # Input checks related to draw functions
         if draw_function:
             # validate function
             validate_draw_function(draw_function, draw_function_kwargs, self.parameters)
-            # function provided but N=1 --> user likely forgot 'N'
-            if N == 1:
-                raise ValueError(
-                    "you specified a draw function but N=1, have you forgotten 'N'?"
-                )
 
         # Copy parameter dictionary --> dict is global
         cp = copy.deepcopy(self.parameters)
@@ -682,7 +670,7 @@ def _mp_sim_single(self, drawn_parameters, time, actual_start_date, method, rtol
         out = self._sim_single(time, actual_start_date, method, rtol, output_timestep, tau)
         return out
 
-    def sim(self, time, warmup=0, N=1, draw_function=None, draw_function_kwargs={}, processes=None, method='RK23', rtol=1e-3, tau=None, output_timestep=1):
+    def sim(self, time, warmup=0, N=1, draw_function=None, draw_function_kwargs={}, processes=None, method='RK23', rtol=1e-4, tau=None, output_timestep=1):
         """
         Simulate a model during a given time period.
         Can optionally perform `N` repeated simulations with sampling of model parameters using a function `draw_function`.
@@ -732,32 +720,15 @@ def sim(self, time, warmup=0, N=1, draw_function=None, draw_function_kwargs={},
         output: xarray.Dataset
             Simulation output
         """
-
-        # Input checks on solver settings
-        if not isinstance(rtol, float):
-            raise TypeError(
-                "relative solver tolerance 'rtol' must be of type float"
-                )
-        if not isinstance(method, str):
-            raise TypeError(
-                "solver method 'method' must be of type string"
-                )
-        if tau != None:
-            if not isinstance(tau, (int,float)):
-                raise TypeError(
-                    "discrete timestep 'tau' must be of type int or float"
-                )
+            
+        # Input checks on solution settings
+        validate_solution_methods_ODE(rtol, method, tau)
         # Input checks on supplied simulation time
         time, actual_start_date = validate_simulation_time(time, warmup)
         # Input checks related to draw functions
         if draw_function:
             # validate function
             validate_draw_function(draw_function, draw_function_kwargs, self.parameters)
-            # function provided but N=1 --> user likely forgot 'N'
-            if N == 1:
-                raise ValueError(
-                    "you specified a draw function but N=1, have you forgotten 'N'?"
-                )
         # provinding 'N' but no draw function: wasteful of resources
         if ((N != 1) & (draw_function==None)):
             raise ValueError('attempting to perform N={0} repeated simulations without using a draw function'.format(N))
diff --git a/src/pySODM/models/validation.py b/src/pySODM/models/validation.py
index ecb487b..2bc9310 100644
--- a/src/pySODM/models/validation.py
+++ b/src/pySODM/models/validation.py
@@ -21,7 +21,7 @@ def validate_simulation_time(time, warmup):
         time = [0-warmup, time]
     elif isinstance(time, list):
         if not len(time) == 2:
-            raise ValueError(f"'Time' must be of format: time=[start, stop]. You have supplied: time={time}.")
+            raise ValueError(f"wrong length of list-like simulation start and stop (length: {len(time)}). correct format: time=[start, stop] (length: 2).")
         else:
             # If they are all int or flat (or commonly occuring np.int64/np.float64)
             if all([isinstance(item, (int,float,np.int32,np.float32,np.int64,np.float64)) for item in time]):
@@ -37,29 +37,87 @@ def validate_simulation_time(time, warmup):
                 actual_start_date = time[0] - timedelta(days=warmup)
                 time = [0, date_to_diff(actual_start_date, time[1])]
             else:
+                types = [type(t) for t in time]
                 raise ValueError(
-                    f"List-like input of simulation start and stop must contain either all int/float or all str/datetime, not a combination of the two "
+                    "simulation start and stop must have the format: time=[start, stop]."
+                    " 'start' and 'stop' must have the same datatype: int/float, str ('yyyy-mm-dd'), or datetime."
+                    f" mixing of types is not allowed. you supplied: {types} "
                     )
-    elif isinstance(time, (str,datetime)):
-        raise TypeError(
-            "You have only provided one date as input 'time', how am I supposed to know when to start/end this simulation?"
-        )
     else:
         raise TypeError(
-                "Input argument 'time' must be a single number (int or float), a list of format: time=[start, stop], a string representing of a timestamp, or a timestamp"
+                "'time' must be 1) a single int/float representing the end of the simulation, 2) a list of format: time=[start, stop]."
             )
 
     if time[1] < time[0]:
         raise ValueError(
-            "Start of simulation is chronologically after end of simulation"
+            "start of simulation is chronologically after end of simulation"
         )
     elif time[0] == time[1]:
-        # TODO: Might be usefull to just return the initial condition in this case?
         raise ValueError(
-            "Start of simulation is the same as the end of simulation"
+            "start of simulation is the same as the end of simulation"
         )
     return time, actual_start_date
 
+def validate_solution_methods_ODE(rtol, method, tau):
+    """
+    Validates the input arguments of the ODE.sim() function
+
+    input
+    -----
+
+    rtol: float
+        Relative solver tolerance
+
+    method: str
+        Solver method: 'RK23', 'RK45', 'DOP853', 'Radau', 'BDF', 'LSODA'
+
+    tau: int/float
+        Discrete integration size of timestep. 
+    """
+
+    if not isinstance(rtol, float):
+        raise TypeError(
+            "relative solver tolerance 'rtol' must be of type float"
+            )
+    if not isinstance(method, str):
+        raise TypeError(
+            "solver method 'method' must be of type string"
+            )
+    if method not in ['RK23', 'RK45', 'DOP853', 'Radau', 'BDF', 'LSODA']:
+        raise ValueError(
+            f"invalid solution method '{method}'. valid methods: 'RK23', 'RK45', 'DOP853', 'Radau', 'BDF', 'LSODA'"
+        )
+    if tau != None:
+        if not isinstance(tau, (int,float)):
+            raise TypeError(
+                "discrete timestep 'tau' must be of type int or float"
+            )
+
+def validate_solution_methods_JumpProcess(method, tau):
+    """
+    Validates the input arguments of the JumpProcess.sim() function
+
+    method: str
+        Solver method: 'SSA' or 'tau_leap'
+
+    tau: int/float
+        If method == 'tau_leap' --> leap size
+    """
+
+    # Input checks on solution method and timestep
+    if not isinstance(method, str):
+        raise TypeError(
+            "solver method 'method' must be of type string"
+            )
+    if method not in ['tau_leap', 'SSA']:
+        raise ValueError(
+            f"invalid solution method '{method}'. valid methods: 'SSA' and 'tau_leap'"
+        )
+    if not isinstance(tau, (int,float)):
+        raise TypeError(
+            "discrete timestep 'tau' must be of type int or float"
+            )
+    
 def validate_draw_function(draw_function, draw_function_kwargs, parameters):
     """Validates the draw function's input and output. For use in the sim() functions of the ODE and JumpProcess classes (base.py).
     
diff --git a/src/pySODM/optimization/objective_functions.py b/src/pySODM/optimization/objective_functions.py
index 3bc8608..78b1709 100644
--- a/src/pySODM/optimization/objective_functions.py
+++ b/src/pySODM/optimization/objective_functions.py
@@ -3,6 +3,7 @@
 import pandas as pd
 import numpy as np
 import xarray as xr
+from datetime import datetime
 from scipy.stats import norm, triang, gamma
 from scipy.special import gammaln
 from pySODM.models.utils import list_to_dict
@@ -268,7 +269,7 @@ class log_posterior_probability():
     # TODO: fully document docstring
     """
     def __init__(self, model, parameter_names, bounds, data, states, log_likelihood_fnc, log_likelihood_fnc_args,
-                 weights=None, log_prior_prob_fnc=None, log_prior_prob_fnc_args=None, initial_states=None, aggregation_function=None, labels=None):
+                 start_sim=None, weights=None, log_prior_prob_fnc=None, log_prior_prob_fnc_args=None, initial_states=None, aggregation_function=None, labels=None):
 
         ############################################################################################################
         ## Validate lengths of data, states, log_likelihood_fnc, log_likelihood_fnc_args, weights, initial states ##
@@ -313,9 +314,12 @@ def __init__(self, model, parameter_names, bounds, data, states, log_likelihood_
 
         # Get additional axes beside time axis in dataset.
         data, self.time_index, self.additional_axes_data = validate_dataset(data)
-        # Extract start- and enddate of simulations
-        self.start_sim = min([df.index.get_level_values(self.time_index).unique().min() for df in data])
-        self.end_sim = max([df.index.get_level_values(self.time_index).unique().max() for df in data])
+
+        ##########################################
+        ## Extract start and end of simulations ##
+        ##########################################
+
+        self.start_sim, self.end_sim = validate_simulation_time_lpp(start_sim, self.time_index, data)
 
         ########################################
         ## Validate the calibrated parameters ##
@@ -518,9 +522,11 @@ def __call__(self, thetas, simulation_kwargs={}):
 def validate_dataset(data):
     """
     Validates a dataset:
-        - Correct datatype?
-        - No Nan's?
-        - Is the index level 'date'/'time present? (obligated)
+        - Does the dataset itself have the right type?
+        - Does it contain Nan's?
+        - Is the index level 'time'/'date' present? (obligated)
+        - Are the indices in 'time' all int/float? Are the indices in 'date' all datetime?
+
     Extracts and returns the additional dimensions in dataset besides the time axis.
 
     Parameters
@@ -535,10 +541,10 @@ def validate_dataset(data):
         List containing the datasets. xarray.DataArray have been converted to pd.DataFrame
 
     additional_axes_data: list
-        Contains the index levels beside 'date'/'time' present in the dataset
+        Contains the index levels beside 'time'/'date' present in the dataset
 
     time_index: str
-        'date': datetime-like time index. 'time': float-like time index.
+        'time': if float-like time index. 'date': if datetime-like time index.
     """
 
     additional_axes_data=[] 
@@ -582,17 +588,97 @@ def validate_dataset(data):
             raise ValueError(
                 f"Index of {idx}th dataset has both 'date' and 'time' as index level (index levels: {df.index.names})."
                 )
-    # Are all time index levels equal to 'date' or 'time'?
+    # Are all time index levels equal to 'date' or 'time' (i.e. no mixing)?
     if all(index == 'date' for index in time_index):
         time_index='date'
     elif all(index == 'time' for index in time_index):
         time_index='time'
     else:
         raise ValueError(
-            "Some datasets have 'time' as time index, other have 'date as time index. pySODM does not allow mixing."
+            "Some datasets have 'time' as temporal index while others have 'date' as temporal index. pySODM does not allow mixing."
         )
+    
+    # Do the types of the time axis make sense?
+    for idx, df in enumerate(data):
+        time_index_vals = df.index.get_level_values(time_index).unique().values
+        # 'time' --> can only be (np.)int/float
+        if time_index=='time':
+            if not all([isinstance(t, (int, float, np.int32, np.int64, np.float32, np.float64)) for t in time_index_vals]):
+                raise ValueError(f"index level 'time' of the {idx}th dataset contains values not of type 'int/float'")
+        # 'date' --> can only be (np.)datetime; no str representations --> gets very messy if we allow this
+        if time_index=='date':
+            if not all([isinstance(t, (datetime, np.datetime64)) for t in time_index_vals]):
+                raise ValueError(f"index level 'date' of the {idx}th dataset contains values not of type 'datetime'")
+
     return data, time_index, additional_axes_data
 
+def validate_simulation_time_lpp(start_sim, time_index, data):
+    """
+    Determines and validates what the start and end of the simulations should be
+
+    Simulation start: user-defined (start_sim is not None) or earliest time/date found in dataasets
+    Simulation end: always latest time/date found in dataasets
+
+    input
+    -----
+
+    start_sim: None (default) or int/float or str/datetime (user-defined)
+        None: start of the simulation is set to the earliest time/date found in dataasets
+        int/float or str/datetime: user-defined start of the simulation
+
+    time_index: str
+        'time': if float-like time index in datasets. 'date': if datetime-like time index in datasets.
+
+    data: list
+        List containing the datasets.
+
+    output
+    ------
+
+    start_sim: float (time_index: 'time') or datetime (time_index: 'date')
+        start of simulations. earliest time/date found in dataasets or user-defined
+
+    end_sim : idem
+        end of simulations. latest time/date found in dataasets
+    """
+
+    # extract startdate
+    if not start_sim:
+        if time_index == 'time':
+            start_sim = float(min([df.index.get_level_values(time_index).unique().min() for df in data]))
+        elif time_index == 'date':
+            start_sim = min([df.index.get_level_values(time_index).unique().min() for df in data]).to_pydatetime() # <class 'pandas._libs.tslibs.timestamps.Timestamp'> --> <class 'datetime.datetime'>
+    else:
+        # assert user input has right datatype
+        assert isinstance(start_sim, (int, float, datetime, str)), "'start_sim' must be of type int, float, str or datetime"
+        # convert to datetime if string
+        if isinstance(start_sim, str):
+            start_sim = datetime.strptime(start_sim,"%Y-%m-%d")
+        # convert to float if int
+        elif isinstance(start_sim, int):
+            start_sim = float(start_sim)
+        else:
+            start_sim = start_sim
+        # only float (time_index: 'time') and datetime (time_index: 'date') remain 
+
+    # extract enddate from datasets
+    if time_index == 'time':
+        end_sim = float(max([df.index.get_level_values(time_index).unique().max() for df in data]))
+    elif time_index == 'date':
+        end_sim = max([df.index.get_level_values(time_index).unique().max() for df in data]).to_pydatetime() 
+
+    # give user pointers on what type to use
+    if time_index == 'time':
+        if not isinstance(start_sim, float):
+            raise TypeError("'start_sim' must be of type int/float")
+    elif time_index == 'date':
+        if not isinstance(start_sim, datetime):
+            raise TypeError("'start_sim' must be of type datetime.datetime or a string representation of a date ('yyyy-mm-dd')")
+
+    # assert startdate before enddate
+    assert start_sim < end_sim, f"'start_sim' ({start_sim}) must be smaller than 'end_sim' ({end_sim})"
+
+    return start_sim, end_sim
 
 def validate_calibrated_parameters(parameters_function, parameters_model):
     """
diff --git a/src/tests/test_JumpProcess.py b/src/tests/test_JumpProcess.py
index 2ab405b..5a52695 100644
--- a/src/tests/test_JumpProcess.py
+++ b/src/tests/test_JumpProcess.py
@@ -1,6 +1,7 @@
 import pytest
-import pandas as pd
 import numpy as np
+import pandas as pd
+from datetime import datetime
 from pySODM.models.base import JumpProcess
 
 #############################
@@ -60,36 +61,39 @@ def test_SIR_time():
     assert S.shape == (51, )
 
     # Do it wrong
+    # dates and times
+    ## Start before end
+    with pytest.raises(ValueError, match="start of simulation is chronologically after end of simulation"):
+        model.sim([20,5])
 
-    # Start same as end
-    with pytest.raises(ValueError, match="Start of simulation is the same as the end of simulation"):
+    ## Start same as end
+    with pytest.raises(ValueError, match="start of simulation is the same as the end of simulation"):
         model.sim([0,0])
 
-    # Start before end
-    with pytest.raises(ValueError, match="Start of simulation is chronologically after end of simulation"):
-        model.sim([20,5])
-
-    # Wrong type
-    with pytest.raises(TypeError, match="Input argument 'time' must be a"):
+    ## Wrong type
+    with pytest.raises(TypeError, match="a single int/float representing the end of the simulation"):
         model.sim(np.zeros(2))    
 
-    # If list: wrong length
-    with pytest.raises(ValueError, match="You have supplied:"):
+    ## If list: wrong length
+    with pytest.raises(ValueError, match="wrong length of list-like simulation start and stop"):
         model.sim([0, 50, 100])    
 
-    # Combination of datetime and int/float
-    with pytest.raises(ValueError, match="List-like input of simulation start and stop"):
-        model.sim([0, pd.to_datetime('2020-03-15')])
+    ## Combination of datetime and int/float
+    with pytest.raises(ValueError, match="simulation start and stop must have the format:"):
+        model.sim([0, datetime(2020,3,15)])
 
-    # Wrong type for input argument `tau`
+    # simulation method
+    ## Wrong type for input argument `tau`
     with pytest.raises(TypeError, match="discrete timestep 'tau' must be of type int or float"):
          model.sim(50, tau='hello')
 
-    # Wrong type for input argument `method`
+    ## Wrong type for input argument `method`
     with pytest.raises(TypeError, match="solver method 'method' must be of type string"):
          model.sim(50, method=3)
 
-test_SIR_time()
+    ## Right type for input argument 'method' but non-existing method
+    with pytest.raises(ValueError, match="invalid solution method 'bliblablu'. valid methods: 'SSA' and 'tau_leap'"):
+        model.sim(50, method='bliblablu')
 
 def test_SIR_date():
     """ Test the use of str/datetime time indexing
@@ -103,7 +107,7 @@ def test_SIR_date():
 
     # Do it right
     output = model.sim(['2020-01-01', '2020-02-20'])
-    output = model.sim([pd.Timestamp('2020-01-01'), pd.Timestamp('2020-02-20')])
+    output = model.sim([datetime(2020,1,1), datetime(2020,2,20)])
 
     # Validate
     assert 'date' in list(output.sizes.keys())
@@ -116,18 +120,26 @@ def test_SIR_date():
     assert S.shape == (51, )
 
     # Do it wrong
-
-    # Start before end
-    with pytest.raises(ValueError, match="Start of simulation is chronologically after end of simulation"):
+    # dates and times
+    ## Start before end
+    with pytest.raises(ValueError, match="start of simulation is chronologically after end of simulation"):
         model.sim(['2020-03-15','2020-03-01'])
 
-    # Combination of str/datetime
-    with pytest.raises(ValueError, match="List-like input of simulation start and stop must contain either"):
-        model.sim([pd.to_datetime('2020-01-01'), '2020-05-01'])
+    ## Start same as end
+    with pytest.raises(ValueError, match="start of simulation is the same as the end of simulation"):
+        model.sim([datetime(2020,1,1),datetime(2020,1,1)])
 
-    # Simulate using a mixture of timestamp and string
-    with pytest.raises(TypeError, match="You have only provided one date as input"):
-        model.sim(pd.to_datetime('2020-01-01'))
+    ## Combination of str/datetime
+    with pytest.raises(ValueError, match="simulation start and stop must have the format:"):
+        model.sim([datetime(2020,1,1), '2020-05-01'])
+
+    ## string not formatted 'yyyy-mm-dd'
+    with pytest.raises(ValueError, match="time data '2020/03/01' does not match format"):
+        model.sim(['2020/03/01', '2020-05-01'])
+
+    ## Only providing one date
+    with pytest.raises(TypeError, match="a single int/float representing the end of the simulation"):
+        model.sim(datetime(2020,1,1))
         
 def test_SSA():
     """ In all testes, the tau-leap is used by default
@@ -140,7 +152,7 @@ def test_SSA():
 
     # Simulate using dates
     output = model.sim(['2020-01-01', '2020-02-20'], method='SSA')
-    output = model.sim([pd.Timestamp('2020-01-01'), pd.Timestamp('2020-02-20')], method='SSA')
+    output = model.sim([datetime(2020,1,1), datetime(2020,2,20)], method='SSA')
 
     # Validate
     assert 'date' in list(output.sizes.keys())
@@ -594,16 +606,12 @@ def draw_function(parameters):
     with pytest.raises(ValueError, match="keys in output dictionary of draw function 'draw_function' must match the keys in input dictionary 'parameters'."):
         model.sim(time, draw_function=draw_function, N=5)
 
-    # correct draw function but user does not provide N
+    # user provides N but no draw function (differs from ODE)
     def draw_function(parameters):
         return parameters
     # simulate model
     time = [0, 10]
     model = SIRstratified(initial_states, parameters, coordinates=coordinates)
-    with pytest.raises(ValueError, match="you specified a draw function but N=1, have you forgotten 'N'"):
-        model.sim(time, draw_function=draw_function)
-
-    # user provides N but no draw function (differs from ODE)
     output = model.sim(time, N=100)
     # assert dimension 'draws' is present in output
     assert 'draws' in list(output.sizes.keys())
diff --git a/src/tests/test_ODE.py b/src/tests/test_ODE.py
index eb27ef9..a6c5594 100644
--- a/src/tests/test_ODE.py
+++ b/src/tests/test_ODE.py
@@ -1,6 +1,7 @@
 import pytest
 import pandas as pd
 import numpy as np
+from datetime import datetime
 from pySODM.models.base import ODE
 
 #############################
@@ -58,28 +59,46 @@ def test_SIR_time():
     assert S.shape == (51, )
 
     # Do it wrong
-    # Start before end
-    with pytest.raises(ValueError, match="Start of simulation is chronologically after end of simulation"):
+    # dates and times
+    ## Start before end
+    with pytest.raises(ValueError, match="start of simulation is chronologically after end of simulation"):
         model.sim([20,5])
 
-    # Start same as end
-    with pytest.raises(ValueError, match="Start of simulation is the same as the end of simulation"):
+    ## Start same as end
+    with pytest.raises(ValueError, match="start of simulation is the same as the end of simulation"):
         model.sim([0,0])
 
-    # Wrong type
-    with pytest.raises(TypeError, match="Input argument 'time' must be a"):
+    ## Wrong type
+    with pytest.raises(TypeError, match="a single int/float representing the end of the simulation"):
         model.sim(np.zeros(2))    
 
-    # If list: wrong length
-    with pytest.raises(ValueError, match="You have supplied:"):
+    ## If list: wrong length
+    with pytest.raises(ValueError, match="wrong length of list-like simulation start and stop"):
         model.sim([0, 50, 100])    
 
-    # Combination of datetime and int/float
-    with pytest.raises(ValueError, match="List-like input of simulation start and stop"):
-        model.sim([0, pd.to_datetime('2020-03-15')])
+    ## Combination of datetime and int/float
+    with pytest.raises(ValueError, match="simulation start and stop must have the format:"):
+        model.sim([0, datetime(2020,3,15)])
+
+    # simulation method
+    ## Wrong type for input argument `rtol`
+    with pytest.raises(TypeError, match="relative solver tolerance 'rtol' must be of type float"):
+         model.sim(50, rtol='hello')
+
+    ## Wrong type for input argument `method`
+    with pytest.raises(TypeError, match="solver method 'method' must be of type string"):
+         model.sim(50, method=3)
+
+    ## Right type for input argument 'method' but non-existing method
+    with pytest.raises(ValueError, match="invalid solution method 'bliblablu'. valid methods:"):
+        model.sim(50, method='bliblablu')
+
+    ## Wrong type for discrete timestep 'tau'
+    with pytest.raises(TypeError, match="discrete timestep 'tau' must be of type int or float"):
+        model.sim(50, tau='bliblablu')
 
 def test_SIR_date():
-    """Test the use of str/datetime time indexing
+    """ Test the use of str/datetime time indexing
     """
 
     # Define parameters and initial states
@@ -90,7 +109,7 @@ def test_SIR_date():
 
     # Do it right
     output = model.sim(['2020-01-01', '2020-02-20'])
-    output = model.sim([pd.Timestamp('2020-01-01'), pd.Timestamp('2020-02-20')])
+    output = model.sim([datetime(2020,1,1), datetime(2020,2,20)])
 
     # Validate
     assert 'date' in list(output.sizes.keys())
@@ -103,19 +122,27 @@ def test_SIR_date():
     assert S.shape == (51, )
 
     # Do it wrong
-
-    # Start before end
-    with pytest.raises(ValueError, match="Start of simulation is chronologically after end of simulation"):
+    # dates and times
+    ## Start before end
+    with pytest.raises(ValueError, match="start of simulation is chronologically after end of simulation"):
         model.sim(['2020-03-15','2020-03-01'])
 
-    # Combination of str/datetime
-    with pytest.raises(ValueError, match="List-like input of simulation start and stop must contain either"):
-        model.sim([pd.to_datetime('2020-01-01'), '2020-05-01'])
+    ## Start same as end
+    with pytest.raises(ValueError, match="start of simulation is the same as the end of simulation"):
+        model.sim([datetime(2020,1,1),datetime(2020,1,1)])
+
+    ## Combination of str/datetime
+    with pytest.raises(ValueError, match="simulation start and stop must have the format:"):
+        model.sim([datetime(2020,1,1), '2020-05-01'])
 
-    # Simulate using a mixture of timestamp and string
-    with pytest.raises(TypeError, match="You have only provided one date as input"):
-        model.sim(pd.to_datetime('2020-01-01'))
+    ## string not formatted 'yyyy-mm-dd'
+    with pytest.raises(ValueError, match="time data '2020/03/01' does not match format"):
+        model.sim(['2020/03/01', '2020-05-01'])
 
+    ## Only providing one date
+    with pytest.raises(TypeError, match="a single int/float representing the end of the simulation"):
+        model.sim(datetime(2020,1,1))
+ 
 def test_SIR_discrete_stepper():
     """ Test the use of the discrete timestepper, `_solve_discrete()`
     """
@@ -571,15 +598,6 @@ def draw_function(parameters):
     with pytest.raises(ValueError, match="keys in output dictionary of draw function 'draw_function' must match the keys in input dictionary 'parameters'."):
         model.sim(time, draw_function=draw_function, N=5)
 
-    # correct draw function but user does not provide N
-    def draw_function(parameters):
-        return parameters
-    # simulate model
-    time = [0, 10]
-    model = SIRstratified(initial_states, parameters, coordinates=coordinates)
-    with pytest.raises(ValueError, match="you specified a draw function but N=1, have you forgotten 'N'"):
-        model.sim(time, draw_function=draw_function)
-
     # user provides N but no draw function
     with pytest.raises(ValueError, match="attempting to perform N=100 repeated simulations without using a draw function"):
         model.sim(time, N=100)
diff --git a/src/tests/test_calibration.py b/src/tests/test_calibration.py
index 3681af2..1ba47fc 100644
--- a/src/tests/test_calibration.py
+++ b/src/tests/test_calibration.py
@@ -1,6 +1,7 @@
 import pytest
-import pandas as pd
 import numpy as np
+import pandas as pd
+from datetime import datetime
 from pySODM.models.base import ODE
 from pySODM.optimization import pso, nelder_mead
 from pySODM.optimization.utils import add_negative_binomial_noise
@@ -20,12 +21,12 @@
 n_datapoints = 20
 alpha = 0.05
 
-###################################################################################
+##############################################################################
 ## Dummy dataset consisting of two dimensions: age groups and spatial units ##
-###################################################################################
+##############################################################################
 
 # Generate a multiindex dataframe with three index levels: time, age groups, spatial componenet
-time = np.linspace(starttime, endtime, num=20)
+time = np.linspace(starttime, endtime, num=n_datapoints)
 age_groups = ['0-20', '20-120']
 spatial_units = [0, 1, 2]
 index = pd.MultiIndex.from_product([time, age_groups, spatial_units], names=['time', 'age_groups', 'spatial_units'])
@@ -70,22 +71,47 @@ def test_weights():
     log_likelihood_fnc_args = [alpha,]
     # Calibated parameters and bounds
     pars = ['beta',]
-    labels = ['$\\beta$',]
     bounds = [(1e-6,1),]
 
     # Correct: list
-    log_posterior_probability(model,pars,bounds,data,states,log_likelihood_fnc,log_likelihood_fnc_args,[1,],labels=labels)
+    log_posterior_probability(model,pars,bounds,data,states,log_likelihood_fnc,log_likelihood_fnc_args,weights=[1,])
     # Correct: np.array
-    log_posterior_probability(model,pars,bounds,data,states,log_likelihood_fnc,log_likelihood_fnc_args,np.array([1,]),labels=labels)
+    log_posterior_probability(model,pars,bounds,data,states,log_likelihood_fnc,log_likelihood_fnc_args,weights=np.array([1,]))
     # Incorrect: list of wrong length
     with pytest.raises(ValueError, match="the extra arguments of the log likelihood function"):
-        log_posterior_probability(model,pars,bounds,data,states,log_likelihood_fnc,log_likelihood_fnc_args,[1,1],labels=labels)
+        log_posterior_probability(model,pars,bounds,data,states,log_likelihood_fnc,log_likelihood_fnc_args,weights=[1,1])
     # Incorrect: numpy array with more than one dimension
     with pytest.raises(TypeError, match="Expected a 1D np.array for input argument"):
-        log_posterior_probability(model,pars,bounds,data,states,log_likelihood_fnc,log_likelihood_fnc_args,np.ones([3,3]),labels=labels)
+        log_posterior_probability(model,pars,bounds,data,states,log_likelihood_fnc,log_likelihood_fnc_args,weights=np.ones([3,3]))
     # Incorrect: datatype string
     with pytest.raises(TypeError, match="Expected a list/1D np.array for input argument"):
-        log_posterior_probability(model,pars,bounds,data,states,log_likelihood_fnc,log_likelihood_fnc_args,'hey',labels=labels)
+        log_posterior_probability(model,pars,bounds,data,states,log_likelihood_fnc,log_likelihood_fnc_args,weights='hey')
+
+def test_start_sim():
+
+    # Define parameters and initial states
+    parameters = {"beta": 0.1, "gamma": 0.2}
+    initial_states = {"S": [1_000_000 - 1], "I": [1], "R": [0]}
+    # Build model
+    model = SIR(initial_states, parameters)
+    # Define dataset
+    data=[df.groupby(by=['time']).sum(),]
+    states = ["I",]
+    log_likelihood_fnc = [ll_negative_binomial,]
+    log_likelihood_fnc_args = [alpha,]
+    # Calibated parameters and bounds
+    pars = ['beta',]
+    bounds = [(1e-6,1),]
+
+    # Correct: int/float (dataset uses 'time' as temporal index)
+    log_posterior_probability(model,pars,bounds,data,states,log_likelihood_fnc,log_likelihood_fnc_args,start_sim=0)
+    # Incorrect: datetime (wrong type)
+    with pytest.raises(TypeError, match="'start_sim' must be of type int/float"):
+        log_posterior_probability(model,pars,bounds,data,states,log_likelihood_fnc,log_likelihood_fnc_args,start_sim=datetime(2020,1,1))  
+    # Incorrect: startdate after enddate
+    with pytest.raises(AssertionError, match="must be smaller than 'end_sim'"):
+        log_posterior_probability(model,pars,bounds,data,states,log_likelihood_fnc,log_likelihood_fnc_args,start_sim=100)  
+
 
 def test_correct_approach_wo_dimension():
 
@@ -106,7 +132,7 @@ def test_correct_approach_wo_dimension():
     bounds = [(1e-6,1),]
     # Setup objective function without priors and with negative weights 
     objective_function = log_posterior_probability(model,pars,bounds,data,states,
-                                                    log_likelihood_fnc,log_likelihood_fnc_args,weights,labels=labels)
+                                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=weights,labels=labels)
 
     # PSO
     theta = pso.optimize(objective_function, kwargs={'simulation_kwargs':{'warmup': warmup}},
@@ -125,7 +151,7 @@ def test_correct_approach_wo_dimension():
         fig_path='sampler_output/'
         identifier = 'username'
         # initialize objective function
-        objective_function = log_posterior_probability(model,pars,bounds,data,states,log_likelihood_fnc,log_likelihood_fnc_args,weights,labels=labels)
+        objective_function = log_posterior_probability(model,pars,bounds,data,states,log_likelihood_fnc,log_likelihood_fnc_args,weights=weights,labels=labels)
         # Perturbate previously obtained estimate
         ndim, nwalkers, pos = perturbate_theta(theta, pert=0.05*np.ones(len(theta)), bounds=objective_function.expanded_bounds, multiplier=multiplier_mcmc)
         # Write some usefull settings to a pickle file (no pd.Timestamps or np.arrays allowed!)
@@ -161,7 +187,7 @@ def break_stuff_wo_dimension():
     # Setup objective function without priors and with negative weights 
     with pytest.raises(Exception, match="is not a valid model parameter!"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)
 
     # Axes in data not present in model
     # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -178,7 +204,7 @@ def break_stuff_wo_dimension():
     # Setup objective function without priors and with negative weights 
     with pytest.raises(Exception, match="Your model has no dimensions."):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)
     
     # Wrong type of dataset
     # ~~~~~~~~~~~~~~~~~~~~~
@@ -196,7 +222,7 @@ def break_stuff_wo_dimension():
     # Setup objective function without priors and with negative weights 
     with pytest.raises(Exception, match="expected pd.Series, pd.DataFrame"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)
 
     # pd.DataFrame with more than one column
     # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -217,7 +243,7 @@ def break_stuff_wo_dimension():
     # Setup objective function without priors and with negative weights 
     with pytest.raises(Exception, match="expected one column."):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)
 
 def break_log_likelihood_functions_wo_dimension():
 
@@ -244,28 +270,28 @@ def break_log_likelihood_functions_wo_dimension():
     # Setup objective function without priors and with negative weights 
     with pytest.raises(ValueError, match="The number of datasets"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)
 
     # wrong type: float
     log_likelihood_fnc = [ll_poisson,]
     log_likelihood_fnc_args = [alpha,]
     with pytest.raises(ValueError, match="dataset has no extra arguments. Expected an empty list"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)
 
     # wrong type: list
     log_likelihood_fnc = [ll_poisson,]
     log_likelihood_fnc_args = [[alpha,],]
     with pytest.raises(ValueError, match="dataset has no extra arguments. Expected an empty list"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)
 
     # wrong type: np.array
     log_likelihood_fnc = [ll_poisson,]
     log_likelihood_fnc_args = [np.array([alpha,]),]
     with pytest.raises(ValueError, match="dataset has no extra arguments. Expected an empty list"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)
 
     # Negative binomial log likelihood
     # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -282,7 +308,7 @@ def break_log_likelihood_functions_wo_dimension():
     log_likelihood_fnc_args = [[alpha,]]
     with pytest.raises(ValueError, match="accepted types are int, float, np.ndarray and pd.Series"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)
 
 def test_xarray_datasets():
     """ Test the use of an xarray.DataArray as dataset
@@ -304,7 +330,7 @@ def test_xarray_datasets():
     bounds = [(1e-6,1),]
     # Setup objective function without priors and with negative weights 
     objective_function = log_posterior_probability(model,pars,bounds,data,states,
-                                                    log_likelihood_fnc,log_likelihood_fnc_args,weights,labels=labels)
+                                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=weights,labels=labels)
 
 class SIR_nd_beta(ODE):
 
@@ -342,7 +368,7 @@ def test_calibration_nd_parameter():
     bounds = [(1e-6,1),]
     # Setup objective function without priors and with negative weights 
     objective_function = log_posterior_probability(model,pars,bounds,data,states,
-                                                    log_likelihood_fnc,log_likelihood_fnc_args,weights,labels=labels)
+                                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=weights,labels=labels)
     # PSO
     theta = pso.optimize(objective_function,
                         swarmsize=10, max_iter=20, processes=1, debug=True)[0]
@@ -393,7 +419,7 @@ def test_correct_approach_with_one_dimension_0():
     data=[df.groupby(by=['time','age_groups']).sum(),]
     # Setup objective function without priors and with negative weights 
     objective_function = log_posterior_probability(model,pars,bounds,data,states,
-                                                    log_likelihood_fnc,log_likelihood_fnc_args,weights,labels=labels)
+                                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=weights,labels=labels)
     # Extract formatted parameter_names, bounds and labels
     labels = objective_function.expanded_labels 
     bounds = objective_function.expanded_bounds
@@ -430,7 +456,7 @@ def break_stuff_with_one_dimension():
     # Setup objective function without priors and with negative weights 
     with pytest.raises(Exception, match="0th dataset coordinate 'spatial_units' is"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)
     
     # Coordinate in dataset not found in the model
     # Setup model
@@ -443,7 +469,7 @@ def break_stuff_with_one_dimension():
     # Setup objective function without priors and with negative weights 
     with pytest.raises(Exception, match="coordinate '0-20' of dimension 'age_groups' in the 0th"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)
 
 def break_log_likelihood_functions_with_one_dimension():
 
@@ -472,27 +498,27 @@ def break_log_likelihood_functions_with_one_dimension():
     log_likelihood_fnc_args = []
     with pytest.raises(ValueError, match="The number of datasets"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)
 
     # wrong type: list of incorrect length
     log_likelihood_fnc = [ll_negative_binomial,]
     log_likelihood_fnc_args = [[alpha,]]
     with pytest.raises(ValueError, match="length of list/1D np.array containing arguments of the log likelihood"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)
 
     # wrong type: np.array of wrong dimensionality
     log_likelihood_fnc = [ll_negative_binomial,]
     log_likelihood_fnc_args = [alpha*np.ones([5,5]), ]
     with pytest.raises(ValueError, match="np.ndarray containing arguments of the log likelihood function"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)    
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)    
 
     # correct type: np.array of right size
     log_likelihood_fnc = [ll_negative_binomial,]
     log_likelihood_fnc_args = [alpha*np.ones([2]), ]
     log_posterior_probability(model,pars,bounds,data,states,
-                                log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)        
+                                log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)        
 
 ####################################
 ## Model with two dimensions ##
@@ -536,7 +562,7 @@ def test_correct_approach_with_two_dimensions():
     data=[df,]
     # Setup objective function without priors and with negative weights 
     objective_function = log_posterior_probability(model,pars,bounds,data,states,
-                                                    log_likelihood_fnc,log_likelihood_fnc_args,weights,labels=labels)
+                                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=weights,labels=labels)
     # Extract formatted parameter_names, bounds and labels
     pars_postprocessing = objective_function.parameters_names_postprocessing
     labels = objective_function.expanded_labels 
@@ -573,14 +599,14 @@ def aggregation_function(output):
         return output.sum(dim='spatial_units')
     # Correct use
     objective_function = log_posterior_probability(model,pars,bounds,data,states,
-                                                    log_likelihood_fnc,log_likelihood_fnc_args,weights,labels=labels,aggregation_function=aggregation_function)
+                                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=weights,labels=labels,aggregation_function=aggregation_function)
     # Correct use
     objective_function = log_posterior_probability(model,pars,bounds,data,states,
-                                                    log_likelihood_fnc,log_likelihood_fnc_args,weights,labels=labels,aggregation_function=[aggregation_function,])
+                                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=weights,labels=labels,aggregation_function=[aggregation_function,])
     # Misuse
     with pytest.raises(ValueError, match="number of aggregation functions must be equal to one or"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,weights,labels=labels,aggregation_function=[aggregation_function,aggregation_function])                                                    
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=weights,labels=labels,aggregation_function=[aggregation_function,aggregation_function])                                                    
 
 def break_log_likelihood_functions_with_two_dimensions():
 
@@ -607,7 +633,7 @@ def break_log_likelihood_functions_with_two_dimensions():
     # Setup objective function without priors and with negative weights 
     with pytest.raises(ValueError, match="Shape of np.array containing arguments of the log likelihood function"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)
     
     # np.array with too many dimensions
     # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -617,7 +643,7 @@ def break_log_likelihood_functions_with_two_dimensions():
     # Setup objective function without priors and with negative weights 
     with pytest.raises(ValueError, match="Shape of np.array containing arguments of the log likelihood function"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)
 
     # np.array with too little dimensions
     # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -627,7 +653,7 @@ def break_log_likelihood_functions_with_two_dimensions():
     # Setup objective function without priors and with negative weights 
     with pytest.raises(ValueError, match="Shape of np.array containing arguments of the log likelihood function"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)
 
     # float
     # ~~~~~
@@ -637,7 +663,7 @@ def break_log_likelihood_functions_with_two_dimensions():
     # Setup objective function without priors and with negative weights 
     with pytest.raises(TypeError, match="accepted types are np.ndarray and pd.Series"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)    
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)    
 
     # float in a list
     # ~~~~~~~~~~~~~~~
@@ -647,7 +673,7 @@ def break_log_likelihood_functions_with_two_dimensions():
     # Setup objective function without priors and with negative weights 
     with pytest.raises(TypeError, match="accepted types are np.ndarray and pd.Series"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)    
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)    
 
     # list
     # ~~~~
@@ -657,7 +683,7 @@ def break_log_likelihood_functions_with_two_dimensions():
     # Setup objective function without priors and with negative weights 
     with pytest.raises(TypeError, match="accepted types are np.ndarray and pd.Series"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)    
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)    
 
     # np.array placed inside too many lists
     # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -667,7 +693,7 @@ def break_log_likelihood_functions_with_two_dimensions():
     # Setup objective function without priors and with negative weights 
     with pytest.raises(TypeError, match="accepted types are np.ndarray and pd.Series"):
         log_posterior_probability(model,pars,bounds,data,states,
-                                    log_likelihood_fnc,log_likelihood_fnc_args,-weights,labels=labels)
+                                    log_likelihood_fnc,log_likelihood_fnc_args,weights=-weights,labels=labels)
 
 ##################################################
 ## Model where states have different dimensions ##
@@ -728,4 +754,4 @@ def test_SIR_SI():
     # Vector population --> calibrate to unstratified data
     data=[df.groupby(by=['time','age_groups']).sum(), df.groupby(by=['time']).sum()]
     # Correct use
-    objective_function = log_posterior_probability(model,pars,bounds,data,states,log_likelihood_fnc,log_likelihood_fnc_args,weights,labels=labels,)
\ No newline at end of file
+    objective_function = log_posterior_probability(model,pars,bounds,data,states,log_likelihood_fnc,log_likelihood_fnc_args,weights=weights,labels=labels,)
\ No newline at end of file
diff --git a/tutorials/influenza_1718/calibration.py b/tutorials/influenza_1718/calibration.py
index 6c5ca15..dff4034 100644
--- a/tutorials/influenza_1718/calibration.py
+++ b/tutorials/influenza_1718/calibration.py
@@ -134,7 +134,6 @@
     # Define dataset
     data=[df_influenza[start_calibration:end_calibration], ]
     states = ["Im_inc",]
-    weights = np.array([1,])
     log_likelihood_fnc = [ll_negative_binomial,]
     log_likelihood_fnc_args = [len(age_groups)*[alpha,],]
     # Calibated parameters and bounds
@@ -142,7 +141,7 @@
     labels = ['$\\beta$', '$f_{ud}$']
     bounds = [(1e-6,0.08), (1e-3,1-1e-3)]
     # Setup objective function (no priors --> uniform priors based on bounds)
-    objective_function = log_posterior_probability(model,pars,bounds,data,states,log_likelihood_fnc,log_likelihood_fnc_args,weights,labels=labels)
+    objective_function = log_posterior_probability(model,pars,bounds,data,states,log_likelihood_fnc,log_likelihood_fnc_args,labels=labels)
     # Extract expanded bounds and labels
     expanded_labels = objective_function.expanded_labels 
     expanded_bounds = objective_function.expanded_bounds