diff --git a/cmake/OpenMCConfig.cmake b/cmake/OpenMCConfig.cmake
index 0bc86fa71c..29a0e4542f 100644
--- a/cmake/OpenMCConfig.cmake
+++ b/cmake/OpenMCConfig.cmake
@@ -1,5 +1,8 @@
 get_filename_component(OpenMC_CMAKE_DIR "${CMAKE_CURRENT_LIST_FILE}" DIRECTORY)
 
+find_package(fmt REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../fmt)
+find_package(gsl-lite REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../gsl-lite)
+find_package(pugixml REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../pugixml)
 find_package(xtl REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtl)
 find_package(xtensor REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtensor)
 
diff --git a/docs/source/usersguide/depletion.rst b/docs/source/usersguide/depletion.rst
index 4a4c7b1f88..f6feb6d245 100644
--- a/docs/source/usersguide/depletion.rst
+++ b/docs/source/usersguide/depletion.rst
@@ -33,12 +33,11 @@ material compositions over time. Each method appears as a different class.
 For example, :class:`openmc.deplete.CECMIntegrator` runs a depletion calculation
 using the CE/CM algorithm (deplete over a timestep using the middle-of-step
 reaction rates). An instance of :class:`openmc.deplete.Operator` is passed to
-one of these functions along with the power level and timesteps::
+one of these functions along with the timesteps and power level::
 
-    power = 1200.0e6
-    days = 24*60*60
-    timesteps = [10.0*days, 10.0*days, 10.0*days]
-    openmc.deplete.CECMIntegrator(op, power, timesteps).integrate()
+    power = 1200.0e6  # watts
+    timesteps = [10.0, 10.0, 10.0]  # days
+    openmc.deplete.CECMIntegrator(op, timesteps, power, timestep_units='d').integrate()
 
 The coupled transport-depletion problem is executed, and once it is done a
 ``depletion_results.h5`` file is written. The results can be analyzed using the
@@ -67,7 +66,7 @@ the energy deposited during a transport calculation will be lower than expected.
 This causes the reaction rates to be over-adjusted to hit the user-specific power,
 or power density, leading to an over-depletion of burnable materials.
 
-There are some remedies. First, the fission Q values can be directly set in a 
+There are some remedies. First, the fission Q values can be directly set in a
 variety of ways. This requires knowing what the total fission energy release should
 be, including indirect components. Some examples are provided below::
 
@@ -99,11 +98,11 @@ Local Spectra and Repeated Materials
 ------------------------------------
 
 It is not uncommon to explicitly create a single burnable material across many locations.
-From a pure transport perspective, there is nothing wrong with creating a single 
+From a pure transport perspective, there is nothing wrong with creating a single
 3.5 wt.% enriched fuel ``fuel_3``, and placing that fuel in every fuel pin in an assembly
 or even full core problem. This certainly expedites the model making process, but can pose
-issues with depletion. 
-Under this setup, :mod:`openmc.deplete` will deplete a single ``fuel_3`` material using 
+issues with depletion.
+Under this setup, :mod:`openmc.deplete` will deplete a single ``fuel_3`` material using
 a single set of reaction rates, and produce a single new composition for the next time
 step. This can be problematic if the same ``fuel_3`` is used in very different regions
 of the problem.
diff --git a/examples/python/pincell_depletion/restart_depletion.py b/examples/python/pincell_depletion/restart_depletion.py
index 013a2469e0..f0387e066f 100644
--- a/examples/python/pincell_depletion/restart_depletion.py
+++ b/examples/python/pincell_depletion/restart_depletion.py
@@ -1,25 +1,7 @@
 import openmc
 import openmc.deplete
-import numpy as np
 import matplotlib.pyplot as plt
 
-###############################################################################
-#                      Simulation Input File Parameters
-###############################################################################
-
-# OpenMC simulation parameters
-batches = 100
-inactive = 10
-particles = 1000
-
-# Depletion simulation parameters
-time_step = 1*24*60*60 # s
-final_time = 5*24*60*60 # s
-time_steps = np.full(final_time // time_step, time_step)
-
-chain_file = './chain_simple.xml'
-power = 174 # W/cm, for 2D simulations only (use W for 3D)
-
 ###############################################################################
 #                      Load previous simulation results
 ###############################################################################
@@ -37,31 +19,34 @@ previous_results = openmc.deplete.ResultsList("depletion_results.h5")
 ###############################################################################
 
 # Instantiate a Settings object, set all runtime parameters
-settings_file = openmc.Settings()
-settings_file.batches = batches
-settings_file.inactive = inactive
-settings_file.particles = particles
+settings = openmc.Settings()
+settings.batches = 100
+settings.inactive = 10
+settings.particles = 10000
 
 # Create an initial uniform spatial source distribution over fissionable zones
 bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1]
 uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
-settings_file.source = openmc.source.Source(space=uniform_dist)
+settings.source = openmc.source.Source(space=uniform_dist)
 
 entropy_mesh = openmc.RegularMesh()
 entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50]
 entropy_mesh.upper_right = [0.39218, 0.39218, 1.e50]
 entropy_mesh.dimension = [10, 10, 1]
-settings_file.entropy_mesh = entropy_mesh
+settings.entropy_mesh = entropy_mesh
 
 ###############################################################################
 #                   Initialize and run depletion calculation
 ###############################################################################
 
-op = openmc.deplete.Operator(geometry, settings_file, chain_file,
-                             previous_results)
+# Create depletion "operator"
+chain_file = './chain_simple.xml'
+op = openmc.deplete.Operator(geometry, settings, chain_file, previous_results)
 
 # Perform simulation using the predictor algorithm
-integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power)
+time_steps = [1.0, 1.0, 1.0, 1.0, 1.0]  # days
+power = 174  # W/cm, for 2D simulations only (use W for 3D)
+integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power, timestep_units='d')
 integrator.integrate()
 
 ###############################################################################
@@ -77,27 +62,28 @@ time, keff = results.get_eigenvalue()
 # Obtain U235 concentration as a function of time
 time, n_U235 = results.get_atoms('1', 'U235')
 
-# Obtain Xe135 absorption as a function of time
-time, Xe_gam = results.get_reaction_rate('1', 'Xe135', '(n,gamma)')
+# Obtain Xe135 capture reaction rate as a function of time
+time, Xe_capture = results.get_reaction_rate('1', 'Xe135', '(n,gamma)')
 
 ###############################################################################
 #                            Generate plots
 ###############################################################################
 
+days = 24*60*60
 plt.figure()
-plt.plot(time/(24*60*60), keff, label="K-effective")
+plt.plot(time/days, keff, label="K-effective")
 plt.xlabel("Time (days)")
 plt.ylabel("Keff")
 plt.show()
 
 plt.figure()
-plt.plot(time/(24*60*60), n_U235, label="U 235")
+plt.plot(time/days, n_U235, label="U 235")
 plt.xlabel("Time (days)")
 plt.ylabel("n U5 (-)")
 plt.show()
 
 plt.figure()
-plt.plot(time/(24*60*60), Xe_gam, label="Xe135 absorption")
+plt.plot(time/days, Xe_capture, label="Xe135 capture")
 plt.xlabel("Time (days)")
 plt.ylabel("RR (-)")
 plt.show()
diff --git a/examples/python/pincell_depletion/run_depletion.py b/examples/python/pincell_depletion/run_depletion.py
index 9933edd48a..ce4b0cf5e6 100644
--- a/examples/python/pincell_depletion/run_depletion.py
+++ b/examples/python/pincell_depletion/run_depletion.py
@@ -1,47 +1,31 @@
+from math import pi
+
 import openmc
 import openmc.deplete
-import numpy as np
 import matplotlib.pyplot as plt
 
-###############################################################################
-#                      Simulation Input File Parameters
-###############################################################################
-
-# OpenMC simulation parameters
-batches = 100
-inactive = 10
-particles = 1000
-
-# Depletion simulation parameters
-time_step = 1*24*60*60 # s
-final_time = 5*24*60*60 # s
-time_steps = np.full(final_time // time_step, time_step)
-chain_file = './chain_simple.xml'
-power = 174 # W/cm, for 2D simulations only (use W for 3D)
-
 ###############################################################################
 #                              Define materials
 ###############################################################################
 
 # Instantiate some Materials and register the appropriate Nuclides
-uo2 = openmc.Material(material_id=1, name='UO2 fuel at 2.4% wt enrichment')
+uo2 = openmc.Material(name='UO2 fuel at 2.4% wt enrichment')
 uo2.set_density('g/cm3', 10.29769)
 uo2.add_element('U', 1., enrichment=2.4)
 uo2.add_element('O', 2.)
-uo2.depletable = True
 
-helium = openmc.Material(material_id=2, name='Helium for gap')
+helium = openmc.Material(name='Helium for gap')
 helium.set_density('g/cm3', 0.001598)
 helium.add_element('He', 2.4044e-4)
 
-zircaloy = openmc.Material(material_id=3, name='Zircaloy 4')
+zircaloy = openmc.Material(name='Zircaloy 4')
 zircaloy.set_density('g/cm3', 6.55)
-zircaloy.add_element('Sn', 0.014  , 'wo')
+zircaloy.add_element('Sn', 0.014, 'wo')
 zircaloy.add_element('Fe', 0.00165, 'wo')
-zircaloy.add_element('Cr', 0.001  , 'wo')
+zircaloy.add_element('Cr', 0.001, 'wo')
 zircaloy.add_element('Zr', 0.98335, 'wo')
 
-borated_water = openmc.Material(material_id=4, name='Borated water')
+borated_water = openmc.Material(name='Borated water')
 borated_water.set_density('g/cm3', 0.740582)
 borated_water.add_element('B', 4.0e-5)
 borated_water.add_element('H', 5.0e-2)
@@ -52,87 +36,62 @@ borated_water.add_s_alpha_beta('c_H_in_H2O')
 #                             Create geometry
 ###############################################################################
 
-# Instantiate ZCylinder surfaces
-fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, r=0.39218, name='Fuel OR')
-clad_ir = openmc.ZCylinder(surface_id=2, x0=0, y0=0, r=0.40005, name='Clad IR')
-clad_or = openmc.ZCylinder(surface_id=3, x0=0, y0=0, r=0.45720, name='Clad OR')
-left = openmc.XPlane(surface_id=4, x0=-0.62992, name='left')
-right = openmc.XPlane(surface_id=5, x0=0.62992, name='right')
-bottom = openmc.YPlane(surface_id=6, y0=-0.62992, name='bottom')
-top = openmc.YPlane(surface_id=7, y0=0.62992, name='top')
+# Define surfaces
+pitch = 1.25984
+fuel_or = openmc.ZCylinder(r=0.39218, name='Fuel OR')
+clad_ir = openmc.ZCylinder(r=0.40005, name='Clad IR')
+clad_or = openmc.ZCylinder(r=0.45720, name='Clad OR')
+box = openmc.model.rectangular_prism(pitch, pitch, boundary_type='reflective')
 
-left.boundary_type = 'reflective'
-right.boundary_type = 'reflective'
-top.boundary_type = 'reflective'
-bottom.boundary_type = 'reflective'
+# Define cells
+fuel = openmc.Cell(fill=uo2, region=-fuel_or)
+gap = openmc.Cell(fill=helium, region=+fuel_or & -clad_ir)
+clad = openmc.Cell(fill=zircaloy, region=+clad_ir & -clad_or)
+water = openmc.Cell(fill=borated_water, region=+clad_or & box)
 
-# Instantiate Cells
-fuel = openmc.Cell(cell_id=1, name='cell 1')
-gap = openmc.Cell(cell_id=2, name='cell 2')
-clad = openmc.Cell(cell_id=3, name='cell 3')
-water = openmc.Cell(cell_id=4, name='cell 4')
-
-# Use surface half-spaces to define regions
-fuel.region = -fuel_or
-gap.region = +fuel_or & -clad_ir
-clad.region = +clad_ir & -clad_or
-water.region = +clad_or & +left & -right & +bottom & -top
-
-# Register Materials with Cells
-fuel.fill = uo2
-gap.fill = helium
-clad.fill = zircaloy
-water.fill = borated_water
-
-# Instantiate Universe
-root = openmc.Universe(universe_id=0, name='root universe')
-
-# Register Cells with Universe
-root.add_cells([fuel, gap, clad, water])
-
-# Instantiate a Geometry, register the root Universe
-geometry = openmc.Geometry(root)
+# Define overall geometry
+geometry = openmc.Geometry([fuel, gap, clad, water])
 
 ###############################################################################
 #                     Set volumes of depletable materials
 ###############################################################################
 
-# Compute cell areas
-area = {}
-area[fuel] = np.pi * fuel_or.coefficients['r'] ** 2
-
-# Set materials volume for depletion. Set to an area for 2D simulations
-uo2.volume = area[fuel]
+# Set material volume for depletion. For 2D simulations, this should be an area.
+uo2.volume = pi * fuel_or.r**2
 
 ###############################################################################
 #                     Transport calculation settings
 ###############################################################################
 
 # Instantiate a Settings object, set all runtime parameters, and export to XML
-settings_file = openmc.Settings()
-settings_file.batches = batches
-settings_file.inactive = inactive
-settings_file.particles = particles
+settings = openmc.Settings()
+settings.batches = 100
+settings.inactive = 10
+settings.particles = 1000
 
 # Create an initial uniform spatial source distribution over fissionable zones
 bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1]
 uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
-settings_file.source = openmc.source.Source(space=uniform_dist)
+settings.source = openmc.source.Source(space=uniform_dist)
 
 entropy_mesh = openmc.RegularMesh()
 entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50]
 entropy_mesh.upper_right = [0.39218, 0.39218, 1.e50]
 entropy_mesh.dimension = [10, 10, 1]
-settings_file.entropy_mesh = entropy_mesh
+settings.entropy_mesh = entropy_mesh
 
 ###############################################################################
 #                   Initialize and run depletion calculation
 ###############################################################################
 
-op = openmc.deplete.Operator(geometry, settings_file, chain_file)
+# Create depletion "operator"
+chain_file = './chain_simple.xml'
+op = openmc.deplete.Operator(geometry, settings, chain_file)
 
 # Perform simulation using the predictor algorithm
-integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power)
+time_steps = [1.0, 1.0, 1.0, 1.0, 1.0]  # days
+power = 174  # W/cm, for 2D simulations only (use W for 3D)
+integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power, timestep_units='d')
 integrator.integrate()
 
 ###############################################################################
@@ -148,27 +107,28 @@ time, keff = results.get_eigenvalue()
 # Obtain U235 concentration as a function of time
 time, n_U235 = results.get_atoms('1', 'U235')
 
-# Obtain Xe135 absorption as a function of time
-time, Xe_gam = results.get_reaction_rate('1', 'Xe135', '(n,gamma)')
+# Obtain Xe135 capture reaction rate as a function of time
+time, Xe_capture = results.get_reaction_rate('1', 'Xe135', '(n,gamma)')
 
 ###############################################################################
 #                            Generate plots
 ###############################################################################
 
+days = 24*60*60
 plt.figure()
-plt.plot(time/(24*60*60), keff, label="K-effective")
+plt.plot(time/days, keff, label="K-effective")
 plt.xlabel("Time (days)")
 plt.ylabel("Keff")
 plt.show()
 
 plt.figure()
-plt.plot(time/(24*60*60), n_U235, label="U 235")
+plt.plot(time/days, n_U235, label="U235")
 plt.xlabel("Time (days)")
 plt.ylabel("n U5 (-)")
 plt.show()
 
 plt.figure()
-plt.plot(time/(24*60*60), Xe_gam, label="Xe135 absorption")
+plt.plot(time/days, Xe_capture, label="Xe135 capture")
 plt.xlabel("Time (days)")
 plt.ylabel("RR (-)")
 plt.show()
diff --git a/include/openmc/particle.h b/include/openmc/particle.h
index 8143ee5a30..2917b45e4a 100644
--- a/include/openmc/particle.h
+++ b/include/openmc/particle.h
@@ -34,12 +34,6 @@ namespace openmc {
 // use to store the bins for delayed group tallies.
 constexpr int MAX_DELAYED_GROUPS {8};
 
-// Maximum number of lost particles
-constexpr int MAX_LOST_PARTICLES {10};
-
-// Maximum number of lost particles, relative to the total number of particles
-constexpr double REL_MAX_LOST_PARTICLES {1.0e-6};
-
 constexpr double CACHE_INVALID {-1.0};
 
 //==============================================================================
diff --git a/include/openmc/settings.h b/include/openmc/settings.h
index 689e222830..494cfd1a20 100644
--- a/include/openmc/settings.h
+++ b/include/openmc/settings.h
@@ -64,10 +64,12 @@ extern std::string path_source_library;   //!< path to the source shared object
 extern std::string path_sourcepoint;      //!< path to a source file
 extern "C" std::string path_statepoint;   //!< path to a statepoint file
 
-extern "C" int32_t n_batches;      //!< number of (inactive+active) batches
-extern "C" int32_t n_inactive;     //!< number of inactive batches
-extern "C" int32_t gen_per_batch;  //!< number of generations per batch
-extern "C" int64_t n_particles;    //!< number of particles per generation
+extern "C" int32_t n_batches;                //!< number of (inactive+active) batches
+extern "C" int32_t n_inactive;               //!< number of inactive batches
+extern "C" int32_t max_lost_particles;     //!< maximum number of lost particles
+extern double rel_max_lost_particles;   //!< maximum number of lost particles, relative to the total number of particles
+extern "C" int32_t gen_per_batch;            //!< number of generations per batch
+extern "C" int64_t n_particles;              //!< number of particles per generation
 
 
 extern int64_t max_particles_in_flight; //!< Max num. event-based particles in flight
diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index 10976dd0ed..31479f9952 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -31,6 +31,10 @@ __all__ = [
     "Integrator", "SIIntegrator", "DepSystemSolver"]
 
 
+_SECONDS_PER_MINUTE = 60
+_SECONDS_PER_HOUR = 60*60
+_SECONDS_PER_DAY = 24*60*60
+
 OperatorResult = namedtuple('OperatorResult', ['k', 'rates'])
 OperatorResult.__doc__ = """\
 Result of applying transport operator
@@ -597,9 +601,11 @@ class Integrator(ABC):
     ----------
     operator : openmc.deplete.TransportOperator
         Operator to perform transport simulations
-    timesteps : iterable of float
-        Array of timesteps in units of [s]. Note that values are not
-        cumulative.
+    timesteps : iterable of float or iterable of tuple
+        Array of timesteps. Note that values are not cumulative. The units are
+        specified by the `timestep_units` argument when `timesteps` is an
+        iterable of float. Alternatively, units can be specified for each step
+        by passing an iterable of (value, unit) tuples.
     power : float or iterable of float, optional
         Power of the reactor in [W]. A single value indicates that
         the power is constant over all timesteps. An iterable
@@ -612,6 +618,11 @@ class Integrator(ABC):
         Power density of the reactor in [W/gHM]. It is multiplied by
         initial heavy metal inventory to get total power if ``power``
         is not speficied.
+    timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
+        Units for values specified in the `timesteps` argument. 's' means
+        seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
+        that the values are given in burnup (MW-d of energy deposited per
+        kilogram of initial heavy metal).
 
     Attributes
     ----------
@@ -625,7 +636,8 @@ class Integrator(ABC):
         Power of the reactor in [W] for each interval in :attr:`timesteps`
     """
 
-    def __init__(self, operator, timesteps, power=None, power_density=None):
+    def __init__(self, operator, timesteps, power=None, power_density=None,
+                 timestep_units='s'):
         # Check number of stages previously used
         if operator.prev_res is not None:
             res = operator.prev_res[-1]
@@ -638,27 +650,59 @@ class Integrator(ABC):
                         self._num_stages))
         self.operator = operator
         self.chain = operator.chain
-        if not isinstance(timesteps, Iterable):
-            self.timesteps = [timesteps]
-        else:
-            self.timesteps = timesteps
+
+        # Determine power and normalize units to W
         if power is None:
             if power_density is None:
                 raise ValueError("Either power or power density must be set")
             if not isinstance(power_density, Iterable):
                 power = power_density * operator.heavy_metal
             else:
-                power = [p * operator.heavy_metal for p in power_density]
-
+                power = [p*operator.heavy_metal for p in power_density]
         if not isinstance(power, Iterable):
             # Ensure that power is single value if that is the case
-            power = [power] * len(self.timesteps)
-        elif len(power) != len(self.timesteps):
-            raise ValueError(
-                "Number of time steps != number of powers. {} vs {}".format(
-                    len(self.timesteps), len(power)))
+            power = [power] * len(timesteps)
 
-        self.power = power
+        if len(power) != len(timesteps):
+            raise ValueError(
+                "Number of time steps ({}) != number of powers ({})".format(
+                    len(timesteps), len(power)))
+
+        # Get list of times / units
+        if isinstance(timesteps[0], Iterable):
+            times, units = zip(*timesteps)
+        else:
+            times = timesteps
+            units = [timestep_units] * len(timesteps)
+
+        # Determine number of seconds for each timestep
+        seconds = []
+        for time, unit, watts in zip(times, units, power):
+            # Make sure values passed make sense
+            check_type('timestep', time, Real)
+            check_greater_than('timestep', time, 0.0, False)
+            check_type('timestep units', unit, str)
+            check_type('power', watts, Real)
+            check_greater_than('power', watts, 0.0, True)
+
+            if unit in ('s', 'sec'):
+                seconds.append(time)
+            elif unit in ('min', 'minute'):
+                seconds.append(time*_SECONDS_PER_MINUTE)
+            elif unit in ('h', 'hr', 'hour'):
+                seconds.append(time*_SECONDS_PER_HOUR)
+            elif unit in ('d', 'day'):
+                seconds.append(time*_SECONDS_PER_DAY)
+            elif unit.lower() == 'mwd/kg':
+                watt_days_per_kg = 1e6*time
+                kilograms = 1e-3*operator.heavy_metal
+                days = watt_days_per_kg * kilograms / watts
+                seconds.append(days*_SECONDS_PER_DAY)
+            else:
+                raise ValueError("Invalid timestep unit '{}'".format(unit))
+
+        self.timesteps = asarray(seconds)
+        self.power = asarray(power)
 
     @abstractmethod
     def __call__(self, conc, rates, dt, power, i):
@@ -772,9 +816,11 @@ class SIIntegrator(Integrator):
     ----------
     operator : openmc.deplete.TransportOperator
         The operator object to simulate on.
-    timesteps : iterable of float
-        Array of timesteps in units of [s]. Note that values are not
-        cumulative.
+    timesteps : iterable of float or iterable of tuple
+        Array of timesteps. Note that values are not cumulative. The units are
+        specified by the `timestep_units` argument when `timesteps` is an
+        iterable of float. Alternatively, units can be specified for each step
+        by passing an iterable of (value, unit) tuples.
     power : float or iterable of float, optional
         Power of the reactor in [W]. A single value indicates that
         the power is constant over all timesteps. An iterable
@@ -787,6 +833,11 @@ class SIIntegrator(Integrator):
         Power density of the reactor in [W/gHM]. It is multiplied by
         initial heavy metal inventory to get total power if ``power``
         is not speficied.
+    timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
+        Units for values specified in the `timesteps` argument. 's' means
+        seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
+        that the values are given in burnup (MW-d of energy deposited per
+        kilogram of initial heavy metal).
     n_steps : int, optional
         Number of stochastic iterations per depletion interval.
         Must be greater than zero. Default : 10
@@ -805,10 +856,10 @@ class SIIntegrator(Integrator):
         Number of stochastic iterations per depletion interval
     """
     def __init__(self, operator, timesteps, power=None, power_density=None,
-                 n_steps=10):
+                 timestep_units='s', n_steps=10):
         check_type("n_steps", n_steps, Integral)
         check_greater_than("n_steps", n_steps, 0)
-        super().__init__(operator, timesteps, power, power_density)
+        super().__init__(operator, timesteps, power, power_density, timestep_units)
         self.n_steps = n_steps
 
     def _get_bos_data_from_operator(self, step_index, step_power, bos_conc):
diff --git a/openmc/deplete/integrators.py b/openmc/deplete/integrators.py
index 67106aa3e1..38a7c33b1f 100644
--- a/openmc/deplete/integrators.py
+++ b/openmc/deplete/integrators.py
@@ -31,9 +31,11 @@ class PredictorIntegrator(Integrator):
     ----------
     operator : openmc.deplete.TransportOperator
         Operator to perform transport simulations
-    timesteps : iterable of float
-        Array of timesteps in units of [s]. Note that values are not
-        cumulative.
+    timesteps : iterable of float or iterable of tuple
+        Array of timesteps. Note that values are not cumulative. The units are
+        specified by the `timestep_units` argument when `timesteps` is an
+        iterable of float. Alternatively, units can be specified for each step
+        by passing an iterable of (value, unit) tuples.
     power : float or iterable of float, optional
         Power of the reactor in [W]. A single value indicates that
         the power is constant over all timesteps. An iterable
@@ -46,6 +48,11 @@ class PredictorIntegrator(Integrator):
         Power density of the reactor in [W/gHM]. It is multiplied by
         initial heavy metal inventory to get total power if ``power``
         is not speficied.
+    timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
+        Units for values specified in the `timesteps` argument. 's' means
+        seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
+        that the values are given in burnup (MW-d of energy deposited per
+        kilogram of initial heavy metal).
 
     Attributes
     ----------
@@ -113,9 +120,11 @@ class CECMIntegrator(Integrator):
     ----------
     operator : openmc.deplete.TransportOperator
         Operator to perform transport simulations
-    timesteps : iterable of float
-        Array of timesteps in units of [s]. Note that values are not
-        cumulative.
+    timesteps : iterable of float or iterable of tuple
+        Array of timesteps. Note that values are not cumulative. The units are
+        specified by the `timestep_units` argument when `timesteps` is an
+        iterable of float. Alternatively, units can be specified for each step
+        by passing an iterable of (value, unit) tuples.
     power : float or iterable of float, optional
         Power of the reactor in [W]. A single value indicates that
         the power is constant over all timesteps. An iterable
@@ -128,6 +137,11 @@ class CECMIntegrator(Integrator):
         Power density of the reactor in [W/gHM]. It is multiplied by
         initial heavy metal inventory to get total power if ``power``
         is not speficied.
+    timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
+        Units for values specified in the `timesteps` argument. 's' means
+        seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
+        that the values are given in burnup (MW-d of energy deposited per
+        kilogram of initial heavy metal).
 
     Attributes
     ----------
@@ -203,9 +217,11 @@ class CF4Integrator(Integrator):
     ----------
     operator : openmc.deplete.TransportOperator
         Operator to perform transport simulations
-    timesteps : iterable of float
-        Array of timesteps in units of [s]. Note that values are not
-        cumulative.
+    timesteps : iterable of float or iterable of tuple
+        Array of timesteps. Note that values are not cumulative. The units are
+        specified by the `timestep_units` argument when `timesteps` is an
+        iterable of float. Alternatively, units can be specified for each step
+        by passing an iterable of (value, unit) tuples.
     power : float or iterable of float, optional
         Power of the reactor in [W]. A single value indicates that
         the power is constant over all timesteps. An iterable
@@ -218,6 +234,11 @@ class CF4Integrator(Integrator):
         Power density of the reactor in [W/gHM]. It is multiplied by
         initial heavy metal inventory to get total power if ``power``
         is not speficied.
+    timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
+        Units for values specified in the `timesteps` argument. 's' means
+        seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
+        that the values are given in burnup (MW-d of energy deposited per
+        kilogram of initial heavy metal).
 
     Attributes
     ----------
@@ -310,9 +331,11 @@ class CELIIntegrator(Integrator):
     ----------
     operator : openmc.deplete.TransportOperator
         Operator to perform transport simulations
-    timesteps : iterable of float
-        Array of timesteps in units of [s]. Note that values are not
-        cumulative.
+    timesteps : iterable of float or iterable of tuple
+        Array of timesteps. Note that values are not cumulative. The units are
+        specified by the `timestep_units` argument when `timesteps` is an
+        iterable of float. Alternatively, units can be specified for each step
+        by passing an iterable of (value, unit) tuples.
     power : float or iterable of float, optional
         Power of the reactor in [W]. A single value indicates that
         the power is constant over all timesteps. An iterable
@@ -325,6 +348,11 @@ class CELIIntegrator(Integrator):
         Power density of the reactor in [W/gHM]. It is multiplied by
         initial heavy metal inventory to get total power if ``power``
         is not speficied.
+    timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
+        Units for values specified in the `timesteps` argument. 's' means
+        seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
+        that the values are given in burnup (MW-d of energy deposited per
+        kilogram of initial heavy metal).
 
     Attributes
     ----------
@@ -404,9 +432,11 @@ class EPCRK4Integrator(Integrator):
     ----------
     operator : openmc.deplete.TransportOperator
         Operator to perform transport simulations
-    timesteps : iterable of float
-        Array of timesteps in units of [s]. Note that values are not
-        cumulative.
+    timesteps : iterable of float or iterable of tuple
+        Array of timesteps. Note that values are not cumulative. The units are
+        specified by the `timestep_units` argument when `timesteps` is an
+        iterable of float. Alternatively, units can be specified for each step
+        by passing an iterable of (value, unit) tuples.
     power : float or iterable of float, optional
         Power of the reactor in [W]. A single value indicates that
         the power is constant over all timesteps. An iterable
@@ -419,6 +449,11 @@ class EPCRK4Integrator(Integrator):
         Power density of the reactor in [W/gHM]. It is multiplied by
         initial heavy metal inventory to get total power if ``power``
         is not speficied.
+    timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
+        Units for values specified in the `timesteps` argument. 's' means
+        seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
+        that the values are given in burnup (MW-d of energy deposited per
+        kilogram of initial heavy metal).
 
     Attributes
     ----------
@@ -518,9 +553,11 @@ class LEQIIntegrator(Integrator):
     ----------
     operator : openmc.deplete.TransportOperator
         Operator to perform transport simulations
-    timesteps : iterable of float
-        Array of timesteps in units of [s]. Note that values are not
-        cumulative.
+    timesteps : iterable of float or iterable of tuple
+        Array of timesteps. Note that values are not cumulative. The units are
+        specified by the `timestep_units` argument when `timesteps` is an
+        iterable of float. Alternatively, units can be specified for each step
+        by passing an iterable of (value, unit) tuples.
     power : float or iterable of float, optional
         Power of the reactor in [W]. A single value indicates that
         the power is constant over all timesteps. An iterable
@@ -533,6 +570,11 @@ class LEQIIntegrator(Integrator):
         Power density of the reactor in [W/gHM]. It is multiplied by
         initial heavy metal inventory to get total power if ``power``
         is not speficied.
+    timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
+        Units for values specified in the `timesteps` argument. 's' means
+        seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
+        that the values are given in burnup (MW-d of energy deposited per
+        kilogram of initial heavy metal).
 
     Attributes
     ----------
@@ -629,9 +671,11 @@ class SICELIIntegrator(SIIntegrator):
     ----------
     operator : openmc.deplete.TransportOperator
         The operator object to simulate on.
-    timesteps : iterable of float
-        Array of timesteps in units of [s]. Note that values are not
-        cumulative.
+    timesteps : iterable of float or iterable of tuple
+        Array of timesteps. Note that values are not cumulative. The units are
+        specified by the `timestep_units` argument when `timesteps` is an
+        iterable of float. Alternatively, units can be specified for each step
+        by passing an iterable of (value, unit) tuples.
     power : float or iterable of float, optional
         Power of the reactor in [W]. A single value indicates that
         the power is constant over all timesteps. An iterable
@@ -644,6 +688,11 @@ class SICELIIntegrator(SIIntegrator):
         Power density of the reactor in [W/gHM]. It is multiplied by
         initial heavy metal inventory to get total power if ``power``
         is not speficied.
+    timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
+        Units for values specified in the `timesteps` argument. 's' means
+        seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
+        that the values are given in burnup (MW-d of energy deposited per
+        kilogram of initial heavy metal).
     n_steps : int, optional
         Number of stochastic iterations per depletion interval.
         Must be greater than zero. Default : 10
@@ -730,9 +779,11 @@ class SILEQIIntegrator(SIIntegrator):
     ----------
     operator : openmc.deplete.TransportOperator
         The operator object to simulate on.
-    timesteps : iterable of float
-        Array of timesteps in units of [s]. Note that values are not
-        cumulative.
+    timesteps : iterable of float or iterable of tuple
+        Array of timesteps. Note that values are not cumulative. The units are
+        specified by the `timestep_units` argument when `timesteps` is an
+        iterable of float. Alternatively, units can be specified for each step
+        by passing an iterable of (value, unit) tuples.
     power : float or iterable of float, optional
         Power of the reactor in [W]. A single value indicates that
         the power is constant over all timesteps. An iterable
@@ -745,6 +796,11 @@ class SILEQIIntegrator(SIIntegrator):
         Power density of the reactor in [W/gHM]. It is multiplied by
         initial heavy metal inventory to get total power if ``power``
         is not speficied.
+    timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
+        Units for values specified in the `timesteps` argument. 's' means
+        seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
+        that the values are given in burnup (MW-d of energy deposited per
+        kilogram of initial heavy metal).
     n_steps : int, optional
         Number of stochastic iterations per depletion interval.
         Must be greater than zero. Default : 10
diff --git a/openmc/lib/settings.py b/openmc/lib/settings.py
index 275c73a934..68d68c29e0 100644
--- a/openmc/lib/settings.py
+++ b/openmc/lib/settings.py
@@ -22,6 +22,8 @@ class _Settings(object):
     entropy_on = _DLLGlobal(c_bool, 'entropy_on')
     generations_per_batch = _DLLGlobal(c_int32, 'gen_per_batch')
     inactive = _DLLGlobal(c_int32, 'n_inactive')
+    max_lost_particles = _DLLGlobal(c_int32, 'max_lost_particles')
+    rel_max_lost_particles = _DLLGlobal(c_double, 'rel_max_lost_particles')
     particles = _DLLGlobal(c_int64, 'n_particles')
     restart_run = _DLLGlobal(c_bool, 'restart_run')
     run_CE = _DLLGlobal(c_bool, 'run_CE')
diff --git a/openmc/settings.py b/openmc/settings.py
index eb2f926f3a..794c97b618 100644
--- a/openmc/settings.py
+++ b/openmc/settings.py
@@ -58,6 +58,10 @@ class Settings(object):
         history-based parallelism.
     generations_per_batch : int
         Number of generations per batch
+    max_lost_particles : int
+        Maximum number of lost particles       
+    rel_max_lost_particles : int
+        Maximum number of lost particles, relative to the total number of particles      
     inactive : int
         Number of inactive batches
     keff_trigger : dict
@@ -176,6 +180,8 @@ class Settings(object):
         self._batches = None
         self._generations_per_batch = None
         self._inactive = None
+        self._max_lost_particles = None
+        self._rel_max_lost_particles = None
         self._particles = None
         self._keff_trigger = None
 
@@ -254,6 +260,14 @@ class Settings(object):
     def inactive(self):
         return self._inactive
 
+    @property
+    def max_lost_particles(self):
+        return self._max_lost_particles        
+
+    @property
+    def rel_max_lost_particles(self):
+        return self._rel_max_lost_particles         
+
     @property
     def particles(self):
         return self._particles
@@ -417,6 +431,19 @@ class Settings(object):
         cv.check_greater_than('inactive batches', inactive, 0, True)
         self._inactive = inactive
 
+    @max_lost_particles.setter
+    def max_lost_particles(self, max_lost_particles):
+        cv.check_type('max_lost_particles', max_lost_particles, Integral)
+        cv.check_greater_than('max_lost_particles', max_lost_particles, 0)
+        self._max_lost_particles = max_lost_particles        
+
+    @rel_max_lost_particles.setter
+    def rel_max_lost_particles(self, rel_max_lost_particles):
+        cv.check_type('rel_max_lost_particles', rel_max_lost_particles, Real)
+        cv.check_greater_than('rel_max_lost_particles', rel_max_lost_particles, 0)
+        cv.check_less_than('rel_max_lost_particles', rel_max_lost_particles, 1)
+        self._rel_max_lost_particles = rel_max_lost_particles          
+
     @particles.setter
     def particles(self, particles):
         cv.check_type('particles', particles, Integral)
@@ -763,6 +790,16 @@ class Settings(object):
             element = ET.SubElement(root, "inactive")
             element.text = str(self._inactive)
 
+    def _create_max_lost_particles_subelement(self, root):
+        if self._max_lost_particles is not None:
+            element = ET.SubElement(root, "max_lost_particles")
+            element.text = str(self._max_lost_particles)    
+
+    def _create_rel_max_lost_particles_subelement(self, root):
+        if self._rel_max_lost_particles is not None:
+            element = ET.SubElement(root, "rel_max_lost_particles")
+            element.text = str(self._rel_max_lost_particles)                     
+
     def _create_particles_subelement(self, root):
         if self._particles is not None:
             element = ET.SubElement(root, "particles")
@@ -1009,6 +1046,8 @@ class Settings(object):
             self._particles_from_xml_element(elem)
             self._batches_from_xml_element(elem)
             self._inactive_from_xml_element(elem)
+            self._max_lost_particles_from_xml_element(elem)
+            self._rel_max_lost_particles_from_xml_element(elem)
             self._generations_per_batch_from_xml_element(elem)
 
     def _run_mode_from_xml_element(self, root):
@@ -1031,6 +1070,16 @@ class Settings(object):
         if text is not None:
             self.inactive = int(text)
 
+    def _max_lost_particles_from_xml_element(self, root):
+        text = get_text(root, 'max_lost_particles')
+        if text is not None:
+            self.max_lost_particles = int(text)     
+
+    def _rel_max_lost_particles_from_xml_element(self, root):
+        text = get_text(root, 'rel_max_lost_particles')
+        if text is not None:
+            self.rel_max_lost_particles = float(text)                       
+
     def _generations_per_batch_from_xml_element(self, root):
         text = get_text(root, 'generations_per_batch')
         if text is not None:
@@ -1270,6 +1319,8 @@ class Settings(object):
         self._create_particles_subelement(root_element)
         self._create_batches_subelement(root_element)
         self._create_inactive_subelement(root_element)
+        self._create_max_lost_particles_subelement(root_element)
+        self._create_rel_max_lost_particles_subelement(root_element)
         self._create_generations_per_batch_subelement(root_element)
         self._create_keff_trigger_subelement(root_element)
         self._create_source_subelement(root_element)
@@ -1340,6 +1391,8 @@ class Settings(object):
         settings._particles_from_xml_element(root)
         settings._batches_from_xml_element(root)
         settings._inactive_from_xml_element(root)
+        settings._max_lost_particles_from_xml_element(root)
+        settings._rel_max_lost_particles_from_xml_element(root)
         settings._generations_per_batch_from_xml_element(root)
         settings._keff_trigger_from_xml_element(root)
         settings._source_from_xml_element(root)
diff --git a/openmc/statepoint.py b/openmc/statepoint.py
index 5802a74cae..84512d70e7 100644
--- a/openmc/statepoint.py
+++ b/openmc/statepoint.py
@@ -310,7 +310,7 @@ class StatePoint(object):
         if self.run_mode == 'eigenvalue':
             return self._f['n_inactive'][()]
         else:
-            return None
+            return None                  
 
     @property
     def n_particles(self):
diff --git a/src/particle.cpp b/src/particle.cpp
index 24bec111c1..d2d8b2e1b3 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -631,8 +631,8 @@ Particle::mark_as_lost(const char* message)
 
   // Abort the simulation if the maximum number of lost particles has been
   // reached
-  if (simulation::n_lost_particles >= MAX_LOST_PARTICLES &&
-      simulation::n_lost_particles >= REL_MAX_LOST_PARTICLES*n) {
+  if (simulation::n_lost_particles >= settings::max_lost_particles &&
+      simulation::n_lost_particles >= settings::rel_max_lost_particles*n) {
     fatal_error("Maximum number of lost particles has been reached.");
   }
 }
diff --git a/src/settings.cpp b/src/settings.cpp
index 6419e46910..614dd8f5b4 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -80,6 +80,8 @@ std::string path_statepoint;
 
 int32_t n_batches;
 int32_t n_inactive {0};
+int32_t max_lost_particles {10};
+double rel_max_lost_particles {1.0e-6};
 int32_t gen_per_batch {1};
 int64_t n_particles {-1};
 
@@ -144,6 +146,16 @@ void get_run_parameters(pugi::xml_node node_base)
   }
   if (!trigger_on) n_max_batches = n_batches;
 
+  // Get max number of lost particles
+  if (check_for_node(node_base, "max_lost_particles")) {
+    max_lost_particles = std::stoi(get_node_value(node_base, "max_lost_particles"));
+  }  
+
+  // Get relative number of lost particles
+  if (check_for_node(node_base, "rel_max_lost_particles")) {
+    rel_max_lost_particles = std::stod(get_node_value(node_base, "rel_max_lost_particles"));
+  }    
+
   // Get number of inactive batches
   if (run_mode == RunMode::EIGENVALUE) {
     if (check_for_node(node_base, "inactive")) {
@@ -344,14 +356,18 @@ void read_settings_xml()
     // Read run parameters
     get_run_parameters(node_mode);
 
-    // Check number of active batches, inactive batches, and particles
+    // Check number of active batches, inactive batches, max lost particles and particles
     if (n_batches <= n_inactive) {
       fatal_error("Number of active batches must be greater than zero.");
     } else if (n_inactive < 0) {
       fatal_error("Number of inactive batches must be non-negative.");
     } else if (n_particles <= 0) {
       fatal_error("Number of particles must be greater than zero.");
-    }
+    } else if (max_lost_particles <= 0) {
+      fatal_error("Number of max lost particles must be greater than zero.");
+    } else if (rel_max_lost_particles <= 0.0 || rel_max_lost_particles >= 1.0) {
+      fatal_error("Relative max lost particles must be between zero and one.");
+    }       
   }
 
   // Copy random number seed if specified
diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py
index 47c9690ae9..6c09b3feac 100644
--- a/tests/unit_tests/test_deplete_integrator.py
+++ b/tests/unit_tests/test_deplete_integrator.py
@@ -7,6 +7,7 @@ will be left unimplemented and testing will be done via regression.
 """
 
 import copy
+from random import uniform
 from unittest.mock import MagicMock
 
 import numpy as np
@@ -15,11 +16,24 @@ import pytest
 
 from openmc.deplete import (
     ReactionRates, Results, ResultsList, comm, OperatorResult,
-    PredictorIntegrator, SICELIIntegrator)
+    PredictorIntegrator, CECMIntegrator, CF4Integrator, CELIIntegrator,
+    EPCRK4Integrator, LEQIIntegrator, SICELIIntegrator, SILEQIIntegrator)
 
 from tests import dummy_operator
 
 
+INTEGRATORS = [
+    PredictorIntegrator,
+    CECMIntegrator,
+    CF4Integrator,
+    CELIIntegrator,
+    EPCRK4Integrator,
+    LEQIIntegrator,
+    SICELIIntegrator,
+    SILEQIIntegrator
+]
+
+
 def test_results_save(run_in_tmpdir):
     """Test data save module"""
 
@@ -105,14 +119,14 @@ def test_results_save(run_in_tmpdir):
     np.testing.assert_array_equal(res[1].time, t2)
 
 
-@pytest.mark.parametrize("timesteps", (1, [1]))
-def test_bad_integrator_inputs(timesteps):
+def test_bad_integrator_inputs():
     """Test failure modes for Integrator inputs"""
 
     op = MagicMock()
     op.prev_res = None
     op.chain = None
     op.heavy_metal = 1.0
+    timesteps = [1]
 
     # No power nor power density given
     with pytest.raises(ValueError, match="Either power or power density"):
@@ -159,3 +173,62 @@ def test_integrator(run_in_tmpdir, scheme):
     dep_time = res.get_depletion_time()
     assert dep_time.shape == (2, )
     assert all(dep_time > 0)
+
+
+@pytest.mark.parametrize("integrator", INTEGRATORS)
+def test_timesteps(integrator):
+    # Crate fake operator
+    op = MagicMock()
+    op.prev_res = None
+    op.chain = None
+
+    # Set heavy metal mass and power randomly
+    op.heavy_metal = uniform(0, 10000)
+    power = uniform(0, 1e6)
+
+    # Reference timesteps in seconds
+    day = 86400.0
+    ref_timesteps = [1*day, 2*day, 5*day, 10*day]
+
+    # Case 1, timesteps in seconds
+    timesteps = ref_timesteps
+    x = integrator(op, timesteps, power, timestep_units='s')
+    assert np.allclose(x.timesteps, ref_timesteps)
+
+    # Case 2, timesteps in minutes
+    minute = 60
+    timesteps = [t / minute for t in ref_timesteps]
+    x = integrator(op, timesteps, power, timestep_units='min')
+    assert np.allclose(x.timesteps, ref_timesteps)
+
+    # Case 3, timesteps in hours
+    hour = 60*60
+    timesteps = [t / hour for t in ref_timesteps]
+    x = integrator(op, timesteps, power, timestep_units='h')
+    assert np.allclose(x.timesteps, ref_timesteps)
+
+    # Case 4, timesteps in days
+    timesteps = [t / day for t in ref_timesteps]
+    x = integrator(op, timesteps, power, timestep_units='d')
+    assert np.allclose(x.timesteps, ref_timesteps)
+
+    # Case 5, timesteps in MWd/kg
+    kilograms = op.heavy_metal / 1000.0
+    days = [t/day for t in ref_timesteps]
+    megawatts = power / 1000000.0
+    burnup = [t * megawatts / kilograms for t in days]
+    x = integrator(op, burnup, power, timestep_units='MWd/kg')
+    assert np.allclose(x.timesteps, ref_timesteps)
+
+    # Case 6, mixed units
+    burnup_per_day = (1e-6*power) / kilograms
+    timesteps = [(burnup_per_day, 'MWd/kg'), (2*day, 's'), (5, 'd'),
+                 (10*burnup_per_day, 'MWd/kg')]
+    x = integrator(op, timesteps, power)
+    assert np.allclose(x.timesteps, ref_timesteps)
+
+    # Bad units should raise an exception
+    with pytest.raises(ValueError, match="unit"):
+        integrator(op, ref_timesteps, power, timestep_units='🐨')
+    with pytest.raises(ValueError, match="unit"):
+        integrator(op, [(800.0, 'gorillas')], power)
diff --git a/tests/unit_tests/test_settings.py b/tests/unit_tests/test_settings.py
index 9c33f6a2d4..c3ff728798 100644
--- a/tests/unit_tests/test_settings.py
+++ b/tests/unit_tests/test_settings.py
@@ -9,6 +9,8 @@ def test_export_to_xml(run_in_tmpdir):
     s.generations_per_batch = 10
     s.inactive = 100
     s.particles = 1000000
+    s.max_lost_particles = 5
+    s.rel_max_lost_particles = 1e-4
     s.keff_trigger = {'type': 'std_dev', 'threshold': 0.001}
     s.energy_mode = 'continuous-energy'
     s.max_order = 5
@@ -62,6 +64,8 @@ def test_export_to_xml(run_in_tmpdir):
     assert s.generations_per_batch == 10
     assert s.inactive == 100
     assert s.particles == 1000000
+    assert s.max_lost_particles == 5
+    assert s.rel_max_lost_particles == 1e-4
     assert s.keff_trigger == {'type': 'std_dev', 'threshold': 0.001}
     assert s.energy_mode == 'continuous-energy'
     assert s.max_order == 5