Steady-state computation in AMICI

This notebook demonstrates how to compute steady states and sensitivities in combination with steady states using AMICI. We explore pre- and post-equilibration scenarios, sensitivity analysis methods, and techniques for handling singular Jacobians.

Steady states play a crucial role in the analysis of dynamical systems, and their computation is a fundamental task in many applications, such as biochemical networks, control systems, and ecological models. A steady state is a condition where the system’s state variables remain constant over time, meaning that their time derivatives are - and remain - zero. In mathematical terms, for a system described by ordinary differential equations (ODEs), a steady state \(\mathbf{x}^*\) satisfies the equation \(\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}^*, \theta) = 0\), where \(\mathbf{f}\) is the right-hand side of the ODEs and \(\theta\) represents model parameters.

AMICI offers several methods for computing steady states and sensitivities at steady state using both forward and adjoint sensitivity analysis.

Steady-state sensitivities are also important in the context of model calibration if the training data contains steady-state measurements and gradient-based optimization methods are used for parameter estimation. Common experimental scenarios include:

1. the system is assumed to start in a steady state, then it is perturbed and enters a dynamic state;

2. the system is assumed to start in a dynamic state, and after some time it reaches a steady state.

In AMICI, these two cases are referred to as pre-equilibration and post-equilibration, respectively. The system can start in a steady state and reach a steady state again after a perturbation, therefore, the two cases are not mutually exclusive.

NOTE: Not every dynamical system needs to run into a steady state. Instead, it may exhibit:

• continuous growth, e.g.,

  \[\dot{x} = x, \quad x_0 = 1\]

• a finite-time blow up, e.g.,

  \[\dot{x} = x^2, \quad x_0 = 1\]

• oscillations, e.g.,

  \[\ddot{x} = -x, \quad x_0 = 1\]

• chaotic behaviour, e.g., the Lorenz attractor

In this notebook, we will discuss how one can perform pre- and post-equilibration using AMICI and compute sensitivities corresponding to the system’s steady state.

The example model

We will use the following model of an enzymatic reaction as example:

The corresponding system of ODEs is

\[\begin{split}\begin{array}{llll} \dot{s} = (\text{synthesis_substrate} + k_{\text{create}}) - k_{\text{bind}}\cdot s \cdot e + k_{\text{unbind}}\cdot c,\\ \dot{e} = - k_{\text{bind}}\cdot s \cdot e + k_{\text{unbind}}\cdot c + k_{\text{convert}}\cdot c,\\ \dot{c} = k_{\text{bind}}\cdot s \cdot e - k_{\text{unbind}}\cdot c - k_{\text{convert}}\cdot c,\\ \dot{p} = k_{\text{convert}}\cdot c - k_{\text{decay}}\cdot p, \end{array}\end{split}\]

where \(s\), \(e\), \(c\) and \(p\) are concentrations of substrate, enzyme, complex and product, respectively.

This model has conserved quantities. Namely, the total amount of enzyme (complex + enzyme) remains constant over time. It can be seen from the ODEs: adding up the second and the third equations gives \(\dot{e} + \dot{c} = 0\). Conserved quantities in the model lead to a singular Jacobian, which makes some computational methods not applicable. For the system above, the Jacobian (\(\mathbf{J} = \frac{\partial \mathbf{f}}{\partial \mathbf{x}}\)) is equal to

\[\begin{split}\begin{bmatrix} - k_{\text{bind}}\cdot e & - k_{\text{bind}}\cdot s & k_{\text{unbind}} & 0\\ - k_{\text{bind}}\cdot e & - k_{\text{bind}}\cdot s & (k_{\text{unbind}} + k_{\text{convert}}) & 0\\ k_{\text{bind}}\cdot e & k_{\text{bind}}\cdot s & -(k_{\text{unbind}} + k_{\text{convert}}) & 0\\ 0 & 0 & k_{\text{convert}} & - k_{\text{decay}}\\ \end{bmatrix}.\end{split}\]

It’s easy to see that the first two columns (or the second and the third row) are linearly dependent and the determinant is equal to zero.

If one is able to identify a conserved quantity, the model dimension can be reduced by excluding one of the states. For example, for the model above one can remove the second ODE and redefine \(e\) as \(\text{enzyme_total} - c\), where \(\text{enzyme_total}\) is a constant specifying the total concentration of the enzyme. In general, conserved quantities or conserved moieties are functions of state variables of the dynamical system that remain constant over time. They can be automatically removed in the SBML import of AMICI, if an environment variable AMICI_EXPERIMENTAL_SBML_NONCONST_CLS is set, using a heuristic-based conserved moieties identification approach presented in De Martino et al. (2014).

In the following, we will consider two versions of this model, with and without eliminating the conserved quantity, to demonstrate the failure of Newton’s method for steady-state computation and of two different approaches for sensitivities computation that are not applicable if the Jacobian is singular.

# Import necessary libraries and define the model
import os
import tempfile
from contextlib import suppress
from pathlib import Path

import amici
import matplotlib.pyplot as plt
import numpy as np
from amici import (
    AMICI_ERROR,
    AMICI_SUCCESS,
    ExpData,
    SensitivityMethod,
    SensitivityOrder,
    SteadyStateComputationMode,
    SteadyStateSensitivityMode,
    SteadyStateStatus,
    import_model_module,
    run_simulation,
    simulation_status_to_str,
)
from amici import (
    MeasurementChannel as MC,
)
from amici.antimony_import import antimony2sbml
from amici.plotting import (
    plot_observable_trajectories,
    plot_state_trajectories,
)

# We encode the model in Antimony format (https://tellurium.readthedocs.io/en/latest/antimony.html),
# which is a human-readable format for biochemical models that can be converted to SBML.
ant_model = """model *model_constant_species()
  const compartment compartment_ = 1;

  // species
  var species substrate in compartment_ = 0;
  var species enzyme in compartment_;
  var species complex in compartment_ = 0;
  species product in compartment_ = 0;

  // dynamic expressions
  var substrate_product := complex + product + substrate;

  // reactions
  creation:  => substrate; compartment_*(synthesis_substrate + k_create);
  binding: substrate + enzyme -> complex; compartment_*(k_bind*substrate*enzyme - k_unbind*complex);
  conversion: complex => enzyme + product; compartment_*k_convert*complex;
  decay: product => ; compartment_*k_decay*product;

  // initial assignments
  enzyme = init_enzyme;
  substrate = init_substrate;

  // parameters
  const init_enzyme = 2;
  const init_substrate = 5;
  const synthesis_substrate = 0;
  const k_bind = 10;
  const k_unbind = 1;
  const k_convert = 10;
  const k_create = 2;
  const k_decay = 1;

  unit volume = 1e-3 litre;
  unit substance = 1e-3 mole;
end"""


# Name of the models that will also be the name of the python module
model_name = "model_constant_species"
model_reduced_name = model_name + "_reduced"

# Directories to which the generated model code is written
temp_dir = Path(tempfile.mkdtemp())
model_output_dir = temp_dir / model_name
model_reduced_output_dir = temp_dir / model_reduced_name

# Import the model
sbml_importer = amici.SbmlImporter(antimony2sbml(ant_model), from_file=False)

# specify observables and constant parameters
constant_parameters = ["synthesis_substrate", "init_enzyme"]
observation_model = [
    MC(id_="observable_product", name="", formula="product"),
    MC(id_="observable_substrate", name="", formula="substrate"),
]
sigmas = {"observable_product": 1.0, "observable_substrate": 1.0}

# import the model without handled the conserved quantity
with suppress(KeyError):
    del os.environ["AMICI_EXPERIMENTAL_SBML_NONCONST_CLS"]
sbml_importer.sbml2amici(
    model_name,
    model_output_dir,
    observation_model=observation_model,
    constant_parameters=constant_parameters,
    compute_conservation_laws=False,
)

# turn on removal of conservation laws and import the reduced model
os.environ["AMICI_EXPERIMENTAL_SBML_NONCONST_CLS"] = "1"
sbml_importer.sbml2amici(
    model_reduced_name,
    model_reduced_output_dir,
    observation_model=observation_model,
    constant_parameters=constant_parameters,
)

<amici.amici.ModelPtr; proxy of <Swig Object of type 'std::unique_ptr< amici::Model > *' at 0x73d676929ef0> >

# import the models and run some test simulations
model_reduced_module = import_model_module(
    model_reduced_name, model_reduced_output_dir
)
model_reduced = model_reduced_module.get_model()

model_module = import_model_module(model_name, model_output_dir)
model = model_module.get_model()

# simulate model with conservation laws
model.set_timepoints(np.linspace(0, 2, 100))
solver = model.create_solver()
rdata = run_simulation(model, solver)
assert rdata.status == AMICI_SUCCESS

# simulate model without conservation laws
model_reduced.set_timepoints(np.linspace(0, 2, 100))
solver_reduced = model_reduced.create_solver()
rdata_reduced = run_simulation(model_reduced, solver_reduced)
assert rdata_reduced.status == AMICI_SUCCESS

# plot trajectories
fig, axes = plt.subplots(2, 2, figsize=(8, 8))
plot_state_trajectories(rdata, model=model, ax=axes[0, 0])
plot_observable_trajectories(rdata, model=model, ax=axes[1, 0])
plot_state_trajectories(rdata_reduced, model=model_reduced, ax=axes[0, 1])
plot_observable_trajectories(rdata_reduced, model=model_reduced, ax=axes[1, 1])
fig.tight_layout()

# the enzyme state was removed from the ODEs of the reduced model
print(model.get_state_ids_solver())
print(model_reduced.get_state_ids_solver())
print(model.get_state_ids())
print(model_reduced.get_state_ids())

('substrate', 'enzyme', 'complex', 'product')
('substrate', 'complex', 'product')
('substrate', 'enzyme', 'complex', 'product')
('substrate', 'enzyme', 'complex', 'product')

Steady-state computation

Two approaches are available in AMICI for steady-state computation:

1. Newton’s method

2. Numerical integration until time derivatives \(\dot{x}\) become sufficiently small. Namely, numerical integration is performed until the condition

   \[\sqrt{\frac 1 n_x \sum_{i=1}^{n_x} (\dot{x}_i w_i)^2} < 1, \quad \text{where} \; w_i = \frac{1}{\text{rtol}*x_i + \text{atol}}\]

   is fulfilled, where “rtol” and “atol” denote relative and absolute tolerances, respectively.

Newton’s method can be significantly faster than numerical integration. However, it will converge only if it was started close enough to the steady state. More about the convergence of the Newton’s method and its implementation in AMICI can be found in Lines et al. (2019). Moreover, Newton’s method is only applicable if the system’s Jacobian is non-singular, and it may converge to a non-physical steady state, i.e., a steady state with negative concentrations (see also Lakrisenko et al. (2024)). Therefore, AMICI users numerical integration by default, which is more robust and can be used for any model.

Many parameters, such as the maximum allowed number of steps for the Newton’s method or numerical integration or values of tolerances, influence steady-state computation.

In the next section, we will try applying Newton’s method to the models with and without conserved quantities, demonstrate how to set relevant parameters and interpret the data returned by simulation.

Inferring the steady state of the system (post-equilibration)

In order to run post-equilibration, one needs to add \(t=\infty\) to the list of timepoints. It can be the only timepoint or one can also set a number of finite timepoints. First, we want to demonstrate that Newton’s method will fail with the unreduced model due to a singular right-hand side Jacobian.

# Call post-equilibration by setting an infinity timepoint
model.set_timepoints([np.inf])

# Set the solver
solver = model.create_solver()
solver.set_newton_max_steps(10)
# maximum number of solver steps for numerical integration
solver.set_max_steps(1000)
model.set_steady_state_computation_mode(
    SteadyStateComputationMode.integrateIfNewtonFails
)
rdata = run_simulation(model, solver)
assert rdata.status == AMICI_SUCCESS
# print out a subset of data returned by model simulation
for key in [
    "ts",
    "x",
    "xdot",
    "posteq_wrms",
    "posteq_t",
    "posteq_numsteps",
    "posteq_numsteps_b",
    "posteq_status",
    "posteq_cpu_time",
    "posteq_cpu_time_b",
]:
    print(f"{key:>16}: {rdata[key]}")

              ts: [inf]
               x: [[0.12222222 1.8        0.2        2.00000314]]
            xdot: [ 0.00000000e+00 -2.22044605e-16  2.22044605e-16 -3.13675608e-06]
     posteq_wrms: 0.7841877852857544
        posteq_t: 14.095659672046486
 posteq_numsteps: [  0 463   0]
posteq_numsteps_b: 0.0
   posteq_status: [<SteadyStateStatus.failed_factorization: -3>, <SteadyStateStatus.success: 1>, <SteadyStateStatus.not_run: 0>]
 posteq_cpu_time: 1.188
posteq_cpu_time_b: 0.0

The x value contains state values corresponding to the time points that were set for the model simulation. In this case, it is the state values at \(t=\infty\), i.e., the steady state. One can also see that xdot values are equal to zero or very small, which means that the steady-state condition (\(\dot{x} \approx 0\)) is fulfilled. The actual weighted root-mean-square of the right-hand side when steady state was reached is reported in posteq_wrms.

The fields posteq_status and posteq_numsteps in rdata tell us how post-equilibration worked. In both lists each entry is, respectively, the status or the number of steps taken for

1. the Newton’s method run;

2. the numerical integration run;

3. the second Newton’s method run, starting from the simulation results.

The status is of type SteadyStateStatus with the following meanings:

• SteadyStateStatus.success: Successful run

• SteadyStateStatus.not_run: Did not run

• SteadyStateStatus.failed: Error: No further specification is given, the error message should give more information.

• SteadyStateStatus.failed_convergence: Error: The method did not converge to a steady state within the maximum number of steps (Newton’s method or simulation).

• SteadyStateStatus.failed_factorization: Error: The Jacobian of the right hand side is singular (only Newton’s method)

• SteadyStateStatus.failed_damping: Error: The damping factor in Newton’s method was reduced until it met the lower bound without success (Newton’s method only)

• SteadyStateStatus.failed_too_long_simulation: Error: The model was simulated past the timepoint t=1e100 without finding a steady state. Therefore, it is likely that the model has no steady state for the given parameter vector.

list(SteadyStateStatus)

[<SteadyStateStatus.failed_too_long_simulation: -5>,
 <SteadyStateStatus.failed_damping: -4>,
 <SteadyStateStatus.failed_factorization: -3>,
 <SteadyStateStatus.failed_convergence: -2>,
 <SteadyStateStatus.failed: -1>,
 <SteadyStateStatus.not_run: 0>,
 <SteadyStateStatus.success: 1>]

In our case, only the second entry of posteq_status contains a positive integer: The first run of Newton’s method failed due to a Jacobian, which could not be factorized, but the second run (simulation) contains the entry 1 (success). The third entry is 0, thus Newton’s method was not launched for the second time. More information can be found inposteq_numsteps: Also here, only the second entry contains a positive integer, which is smaller than the allowed maximum number of steps (1000). Hence, steady state was reached via simulation, which corresponds to the simulated time written to posteq_time.

We can simulate the system until \(t=14\) and see that, indeed, the system converges to the steady-state value, which is equal to rdata['x'][0].

# steady state value found by pre-equilibration
steady_state = rdata["x"][0]

timepoints = np.linspace(0, 14, 200)
model.set_timepoints(timepoints)
rdata = run_simulation(model, solver)
plot_state_trajectories(rdata, model=model)

for stst_value in steady_state:
    plt.axhline(y=stst_value, color="gray", linestyle="--", linewidth=1)

We want to demonstrate a complete failure during the steady-state computation by reducing the number of integration steps to a lower value:

# reduce maxsteps for integration
model.set_timepoints([np.inf])
solver.set_max_steps(100)
rdata = run_simulation(model, solver)
assert rdata.status == AMICI_ERROR
print("Simulation status:", simulation_status_to_str(rdata["status"]))
print("Status of post-equilibration:", rdata["posteq_status"])
print(
    "Number of steps employed in post-equilibration:", rdata["posteq_numsteps"]
)

2025-11-10 08:36:13.150 - amici.swig_wrappers - DEBUG - [EQUILIBRATION_FAILURE] AMICI equilibration exceeded maximum number of integration steps at t=0.0210298.
2025-11-10 08:36:13.150 - amici.swig_wrappers - ERROR - [OTHER] AMICI simulation failed: Steady state computation failed. First run of Newton solver failed: RHS could not be factorized. Simulation to steady state failed: No convergence was achieved. Second run of Newton solver failed: RHS could not be factorized.

Simulation status: AMICI_ERROR
Status of post-equilibration: [<SteadyStateStatus.failed_factorization: -3>, <SteadyStateStatus.failed_convergence: -2>, <SteadyStateStatus.failed_factorization: -3>]
Number of steps employed in post-equilibration: [  0 100   0]

However, the same logic works if we use the reduced model. For sufficiently many Newton steps, post-equilibration is achieved by Newton’s method in the first run. In this specific example, the steady state is found within one step.

model_reduced.set_timepoints([np.inf])
model_reduced.set_steady_state_computation_mode(
    SteadyStateComputationMode.integrateIfNewtonFails
)

# set the solver
solver_reduced = model_reduced.create_solver()
solver_reduced.set_newton_max_steps(10)
solver_reduced.set_max_steps(100)
rdata_reduced = run_simulation(model_reduced, solver_reduced)
assert rdata_reduced.status == AMICI_SUCCESS
assert rdata_reduced.posteq_status[0] == SteadyStateStatus.success
print(
    "Simulation status:",
    simulation_status_to_str(rdata_reduced["status"]),
)
print("Steady-state value:", rdata_reduced["x"][0])
print("Status of post-equilibration:", rdata_reduced["posteq_status"])
print(
    "Number of steps employed in post-equilibration:",
    rdata_reduced["posteq_numsteps"],
)

Simulation status: AMICI_SUCCESS
Steady-state value: [0.12222222 1.8        0.2        2.        ]
Status of post-equilibration: [<SteadyStateStatus.success: 1>, <SteadyStateStatus.not_run: 0>, <SteadyStateStatus.not_run: 0>]
Number of steps employed in post-equilibration: [2 0 0]

Pre-equilibrating the model

Sometimes, we want to launch a solver run from a steady state which was inferred numerically, i.e., the system requires pre-equilibration. In order to do this with AMICI, we need to pass an ExpData object, which contains fixed parameter values for the actual simulation and for pre-equilibration of the model.

# create edata, with 3 timepoints and 2 observables:
edata = ExpData(2, 0, 0, np.array([0.0, 0.1, 1.0]))
edata.set_observed_data([1.8] * 6)
edata.fixed_parameters = np.array([3.0, 5.0])
# set parameters for pre-equilibration
edata.fixed_parameters_pre_equilibration = np.array([0.0, 2.0])
edata.reinitialize_fixed_parameter_initial_states = True

# create the solver object and run the simulation
solver_reduced = model_reduced.create_solver()
solver_reduced.set_newton_max_steps(10)
rdata_reduced = run_simulation(model_reduced, solver_reduced, edata)
assert rdata_reduced.status == AMICI_SUCCESS
print(
    "Simulation status:",
    simulation_status_to_str(rdata_reduced["status"]),
)
print("Status of pre-equilibration:", rdata_reduced["preeq_status"])
print(
    "Number of steps employed in pre-equilibration:",
    rdata_reduced["preeq_numsteps"],
)
print("Steady state found during pre-equilibration:", rdata_reduced["x_ss"])
print("Initial state used for simulation:", rdata_reduced["x0"])


plot_state_trajectories(rdata_reduced, model=model_reduced)
plot_observable_trajectories(rdata_reduced, model=model_reduced)

Simulation status: AMICI_SUCCESS
Status of pre-equilibration: [<SteadyStateStatus.success: 1>, <SteadyStateStatus.not_run: 0>, <SteadyStateStatus.not_run: 0>]
Number of steps employed in pre-equilibration: [2 0 0]
Steady state found during pre-equilibration: [0.12222222 1.8        0.2        2.        ]
Initial state used for simulation: [0.12222222 5.         0.2        2.        ]

Pre-equilibration was successful and required 2 steps of Newton’s method. The steady state found during pre-equilibration is used as the initial state for the subsequent simulation, except the second entry, corresponding to the enzyme concentration. The reason is that the initial value for the enzyme concentration is defined by the second fixed parameter (init_enzyme). This parameter value is re-initialised after pre-equilibration and is set to 5.

We can also combine pre- and post-equilibration:

# Change the last timepoint to an infinity timepoint.
edata.set_timepoints(np.array([0.0, 0.1, float("inf")]))

# run the simulation
rdata_reduced = run_simulation(model_reduced, solver_reduced, edata)
assert rdata_reduced.status == AMICI_SUCCESS

Sensitivities computation at steady state

In addition to model simulation, AMICI supports forward sensitivity analysis and adjoint sensitivity analysis for likelihood-based output functions, which is crucial for parameter estimation. Both sensitivity analysis approaches can be adapted to exploit steady-state constraints and to compute sensitivities by solving a linear system. In this section, we consider these approaches, their applicability and available settings.

Post-equilibration with sensitivities

Post-equilibration with forward sensitivities

If forward sensitivity analysis is used, then state sensitivities at the timepoint np.inf will be computed. This can be done in (currently) two different ways:

1. One can integrate the state variables with state sensitivities until the norm of the right-hand side becomes small. This approach is numerically more stable, but the computation time for large models may be substantial.

2. The more efficient approach is to take into account the steady-state condition (\(\dot{x} = 0\)), which allows to simplify the ODE system for state sensitivities into the linear system of equations:

   \[ 0 = \left.\dot{\mathbf{s}}^x\right\vert_{\mathbf{x} = \mathbf{x}^*({\theta})} = \mathbf{J}(\mathbf{x}^*({\theta}), {\theta}) \left.\mathbf{s}^x\right\vert_{\mathbf{x} = \mathbf{x}^*({\theta})} + \left.\frac{\partial \mathbf{f}}{\partial \theta}\right\vert_{\mathbf{x} = \mathbf{x}^*({\theta})} \qquad \Rightarrow \qquad \mathbf{J}(\mathbf{x}^*({\theta}), {\theta}) \left.\mathbf{s}^x\right\vert_{\mathbf{x} = \mathbf{x}^*({\theta})} = - \left.\frac{\partial \mathbf{f}}{\partial \theta}\right\vert_{\mathbf{x} = \mathbf{x}^*({\theta})},\]

   where \(\mathbf{s}^x\) are the state sensitivities, \(\mathbf{f}\) is the right-hand side of the model ODEs, \(\theta\) the model parameters, and \(\mathbf{J}(\mathbf{x}^*({\theta}), {\theta}) = \left.\nabla_x \mathbf{f}\right\vert_{\mathbf{x} = \mathbf{x}^*({\theta})}\) is the systems’s Jacobian at steady state. However, this method is only applicable if the Jacobian is not (close to) singular.

This approach will always be chosen by AMICI, if the steady-state computation mode is set to SteadyStateSensitivityMode.newtonOnly (via Model.setSteadyStateSensitivityMode()). Furthermore, it will also be chosen if the steady state was found by Newton’s method, as in this case, the Jacobian is at least not singular (but may still be poorly conditioned).

Side remark: A possible third way may consist in a (relaxed) Richardson iteration type approach, which interprets the entries of the right hand side \(f\) as residuals and minimizes the squared residuals \(\Vert f \Vert^2\) by a Levenberg-Marquardt-type algorithm. This approach would also work for poorly conditioned (and even for singular Jacobians if additional constraints are implemented as Lagrange multipliers) while being faster than a long forward simulation.

We want to demonstrate both possibilities for steady-state sensitivities computation, and the failure of their computation if the Jacobian is singular and the newtonOnly setting was used.

# Call simulation with singular Jacobian and `integrateIfNewtonFails` mode
model.set_timepoints([np.inf])
model.set_steady_state_sensitivity_mode(
    SteadyStateSensitivityMode.integrateIfNewtonFails
)
solver = model.create_solver()
solver.set_newton_max_steps(10)
solver.set_sensitivity_method(SensitivityMethod.forward)
solver.set_sensitivity_order(SensitivityOrder.first)
solver.set_max_steps(10000)
rdata = run_simulation(model, solver)

np.set_printoptions(threshold=20)
print("Simulation status:", simulation_status_to_str(rdata["status"]))
print("Status of post-equilibration:", rdata["posteq_status"])
print(
    "Number of steps employed in post-equilibration:", rdata["posteq_numsteps"]
)
print("Computed state sensitivities:")
print(rdata["sx"][0, :, :])

Simulation status: AMICI_SUCCESS
Status of post-equilibration: [<SteadyStateStatus.failed_factorization: -3>, <SteadyStateStatus.success: 1>, <SteadyStateStatus.not_run: 0>]
Number of steps employed in post-equilibration: [   0 1346    0]
Computed state sensitivities:
[[ 6.79012346e-02 -1.00000000e-01  1.00000000e-01  1.00000000e+00]
 [-1.22222222e-02 -3.50756789e-17  4.45553752e-18 -4.96719843e-16]
 [ 1.11111111e-02 -6.56873594e-18  4.89709017e-18  4.77663120e-16]
 [-2.46913580e-03  2.00000000e-02 -2.00000000e-02 -1.28543060e-15]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00 -2.00000000e+00]
 [ 6.17862181e-19 -9.09942485e-18 -1.34718585e-34  1.67864932e-14]]

As can be seen from the posteq_status, Newton’s method didn’t work, but the steady state and corresponding sensitivities were computed by numerical integration. The shape of rdata['sx'] is (number of timepoints) x (number of parameters) x (number of states).

# Call simulation with singular Jacobian and newtonOnly mode (will fail)
model.set_timepoints([np.inf])
model.set_steady_state_sensitivity_mode(SteadyStateSensitivityMode.newtonOnly)
solver = model.create_solver()
solver.set_newton_max_steps(10)
solver.set_sensitivity_method(SensitivityMethod.forward)
solver.set_sensitivity_order(SensitivityOrder.first)
solver.set_max_steps(10000)
rdata = run_simulation(model, solver)
assert rdata.status == AMICI_ERROR

print("Simulation status:", simulation_status_to_str(rdata["status"]))
print("Steady state:", rdata["x"])
print("Status of post-equilibration:", rdata["posteq_status"])
print(
    "Number of steps employed in post-equilibration:", rdata["posteq_numsteps"]
)
print("Computed state sensitivities:")
print(rdata["sx"][0, :, :])

2025-11-10 08:36:13.709 - amici.swig_wrappers - ERROR - [OTHER] AMICI simulation failed: Steady state sensitivity computation failed due to unsuccessful factorization of RHS Jacobian

Simulation status: AMICI_ERROR
Steady state: [[0.12222222 1.8        0.2        2.00000314]]
Status of post-equilibration: [<SteadyStateStatus.failed_factorization: -3>, <SteadyStateStatus.success: 1>, <SteadyStateStatus.not_run: 0>]
Number of steps employed in post-equilibration: [  0 463   0]
Computed state sensitivities:
[[0. 0. 0. 0.]
 [0. 0. 0. 0.]
 [0. 0. 0. 0.]
 [0. 0. 0. 0.]
 [0. 0. 0. 0.]
 [1. 0. 0. 0.]]

For the model with singular Jacobian steady-state computation worked by numerical integration. However, for sensitivities computation SteadyStateSensitivityMode.newtonOnly approach was chosen which fails with singular Jacobian and the sensitivities could not be computed. We can try the same with the reduced model:

# Try `newtonOnly` option with reduced model
model_reduced.set_timepoints([np.inf])
model_reduced.set_steady_state_sensitivity_mode(
    SteadyStateSensitivityMode.newtonOnly
)
solver_reduced = model_reduced.create_solver()
solver_reduced.set_newton_max_steps(10)
solver_reduced.set_sensitivity_method(SensitivityMethod.forward)
solver_reduced.set_sensitivity_order(SensitivityOrder.first)
solver_reduced.set_max_steps(1000)
rdata_reduced = run_simulation(model_reduced, solver_reduced)
assert rdata_reduced.status == AMICI_SUCCESS
assert rdata_reduced.posteq_status == [
    SteadyStateStatus.success,
    SteadyStateStatus.not_run,
    SteadyStateStatus.not_run,
]

print(
    "Simulation status:",
    simulation_status_to_str(rdata_reduced["status"]),
)
print("Steady state:", rdata_reduced["x"])
print("Status of post-equilibration:", rdata_reduced["posteq_status"])
print(
    "Number of steps employed in post-equilibration:",
    rdata_reduced["posteq_numsteps"],
)
print("Computed state sensitivities:")
print(rdata_reduced["sx"][0, :, :])

Simulation status: AMICI_SUCCESS
Steady state: [[0.12222222 1.8        0.2        2.        ]]
Status of post-equilibration: [<SteadyStateStatus.success: 1>, <SteadyStateStatus.not_run: 0>, <SteadyStateStatus.not_run: 0>]
Number of steps employed in post-equilibration: [2 0 0]
Computed state sensitivities:
[[ 0.06790123 -0.1         0.1         1.        ]
 [-0.01222222  0.          0.          0.        ]
 [ 0.01111111  0.          0.          0.        ]
 [-0.00246914  0.02       -0.02        0.        ]
 [ 0.          0.          0.         -2.        ]
 [ 0.          0.          0.          0.        ]]

In this case, both steady state and the corresponding sensitivities could be computed by solving the linear systems \(\mathbf{s}^x = - \frac{\partial \mathbf{f}}{\partial \theta}\) for state sensitivities (\(\mathbf{s}^x\)).

Post-equilibration with adjoint sensitivities

Post-equilibration also works with adjoint sensitivity analysis and, similar to forward sensitivities, two approaches are available for sensitivities computation at steady state:

1. Standard adjoint sensitivity analysis can be performed to compute the objective function gradient. This approach requires backward numerical integration of the adjoint state ODE form the timepoint of the steady-state measurements \(t=t''\) (practically, \(t''\) = posteq_t, at this time point steady-state condition was satisfied) until the last finite timepoint (\(t_{n_t}\)). Ultimately, one needs to compute

   \[\int_{t_{n_t}}^{t''}\mathbf{p}(t, {\theta})^T\left.\frac{\partial \mathbf{f}}{\partial \theta_k}\right|_{\textbf{x}(t, {\theta}),{\theta}}dt,\]

   where \(\mathbf{p}\) is the adjoint state, \(\theta\) the model parameters.

2. Alternatively, one can exploit the steady-state condition and speed up the computation. It can be shown that the integral

   \[\int_{t_{n_t}}^{t''}\mathbf{p}(t, {\theta})^T\left.\frac{\partial \mathbf{f}}{\partial \theta_k}\right|_{\textbf{x}(t, {\theta}),{\theta}}dt\]

   can be reduced to a matrix-vector product

   \[\int^{t''}_{t'} \mathbf{p}(t, {\theta}, \mathbf{u})^T dt \cdot \left. \frac{\partial \mathbf{f}} {\partial \theta_k} \right\vert_{\mathbf{x} = \mathbf{x}^*({\theta})}dt = \mathbf{p}_{\text{integral}} \cdot \left. \frac{\partial \mathbf{f}} {\partial \theta_k} \right\vert_{\mathbf{x} = \mathbf{x}^*({\theta})},\]

   where \(t'\) is a time point such that \(t_{n_t} \ll t' \ll t''\) and for \(t \geq t'\) the system is at steady state. To obtain the integral one can solve the linear system of equations

   \[\mathbf{J}(\mathbf{x}^*({\theta}), {\theta})^T \mathbf{p}_{\text{integral}} = -\mathbf{p}(t'', {\theta}).\]

   The detailed explanation of this approach and the derivation of the last equation can be found in Lakrisenko et al. (2023).

   However, this solution is given in terms of a linear system of equations defined by the transposed Jacobian of the right-hand side. Hence, if the (transposed) Jacobian is singular, it is not applicable. In this case, standard integration must be carried out.

NOTE: adjoint sensitivity analysis yields less information than forward sensitivity analysis, since state sensitivities are not computed. However, it has a better scaling behavior towards large model size.

# Call post-equilibration and sensitivities computation using adjoint sensitivity analysis
# by setting an infinity timepoint
# and creatíng an edata object, which is needed for adjoint computation
edata = ExpData(2, 0, 0, np.array([float("inf")]))
edata.set_observed_data([1.8] * 2)
edata.fixed_parameters = np.array([3.0, 5.0])

model_reduced.set_steady_state_sensitivity_mode(
    SteadyStateSensitivityMode.newtonOnly
)
solver_reduced = model_reduced.create_solver()
solver_reduced.set_newton_max_steps(10)
solver_reduced.set_sensitivity_method(SensitivityMethod.adjoint)
solver_reduced.set_sensitivity_order(SensitivityOrder.first)
solver_reduced.set_max_steps(1000)
rdata_reduced = run_simulation(model_reduced, solver_reduced, edata)

print(
    "Simulation status:",
    simulation_status_to_str(rdata_reduced["status"]),
)
print("Status of post-equilibration:", rdata_reduced["posteq_status"])
print(
    "Number of steps employed in post-equilibration:",
    rdata_reduced["posteq_numsteps"],
)
print(
    "Number of backward steps employed in post-equilibration:",
    rdata_reduced["posteq_numsteps_b"],
)
print("Computed gradient:", rdata_reduced["sllh"])

Simulation status: AMICI_SUCCESS
Status of post-equilibration: [<SteadyStateStatus.success: 1>, <SteadyStateStatus.not_run: 0>, <SteadyStateStatus.not_run: 0>]
Number of steps employed in post-equilibration: [2 0 0]
Number of backward steps employed in post-equilibration: 0.0
Computed gradient: [-3.15443073e+00 -2.05061728e-02  1.86419753e-02 -4.14266118e-03
  1.60000000e+01  0.00000000e+00]

If we carry out the same computation with a system that has a singular Jacobian, then posteq_numsteps_b will not be 0 anymore (which indicates that the linear system solve was used to compute backward post-equilibration). Now, integration is carried out and hence posteq_numsteps_b > 0:

# Call adjoint post-equilibration with model with singular Jacobian
model.set_steady_state_sensitivity_mode(SteadyStateSensitivityMode.newtonOnly)
solver = model.create_solver()
solver.set_newton_max_steps(10)
solver.set_max_steps_backward_problem(10000)
solver.set_sensitivity_method(SensitivityMethod.adjoint)
solver.set_sensitivity_order(SensitivityOrder.first)
rdata = run_simulation(model, solver, edata)

print("Simulation status:", simulation_status_to_str(rdata["status"]))
print("Status of post-equilibration:", rdata["posteq_status"])
print(
    "Number of steps employed in post-equilibration:", rdata["posteq_numsteps"]
)
print(
    "Number of backward steps employed in post-equilibration:",
    rdata["posteq_numsteps_b"],
)
print("Computed gradient:", rdata["sllh"])

2025-11-10 08:36:13.733 - amici.swig_wrappers - ERROR - [OTHER] AMICI simulation failed: Steady state backward computation failed: Linear system could not be solved (possibly due to singular Jacobian), and numerical integration did not equilibrate within maxsteps

Simulation status: AMICI_ERROR
Status of post-equilibration: [<SteadyStateStatus.failed_factorization: -3>, <SteadyStateStatus.success: 1>, <SteadyStateStatus.not_run: 0>]
Number of steps employed in post-equilibration: [  0 408   0]
Number of backward steps employed in post-equilibration: 0.0
Computed gradient: [nan nan nan nan nan nan]

Pre-equilibration with sensitivities

Beyond the need for an ExpData object, the steady state solver logic in pre-equilibration is the same as in post-equilibration, also if sensitivities are requested. The computation will fail for singular Jacobians, if SteadyStateSensitivityMode is set to newtonOnly, or if not enough steps can be taken. However, if forward simulation with steady state sensitivities is allowed, or if the Jacobian is not singular, it will work.

Pre-equilibration with forward sensitivities

# No post-equilibration this time.
# create edata, with 3 timepoints and 2 observables:
edata = ExpData(2, 0, 0, np.array([0.0, 0.1, 1.0]))
edata.set_observed_data([1.8] * 6)
edata.fixed_parameters = np.array([3.0, 5.0])
# set parameters for pre-equilibration
edata.fixed_parameters_pre_equilibration = np.array([0.0, 2.0])
edata.reinitialize_fixed_parameter_initial_states = True

# create the solver object and run the simulation, singular Jacobian, enforce Newton solver for sensitivities
model.set_steady_state_sensitivity_mode(SteadyStateSensitivityMode.newtonOnly)
solver = model.create_solver()
solver.set_newton_max_steps(10)
solver.set_sensitivity_method(SensitivityMethod.forward)
solver.set_sensitivity_order(SensitivityOrder.first)
rdata = run_simulation(model, solver, edata)

assert rdata.status == AMICI_ERROR
for key, value in rdata.items():
    if key[0:6] == "preeq_":
        print(f"{key:20s}:", value)

2025-11-10 08:36:13.743 - amici.swig_wrappers - ERROR - [OTHER] AMICI simulation failed: Steady state sensitivity computation failed due to unsuccessful factorization of RHS Jacobian

preeq_wrms          : 0.7841877852857544
preeq_t             : 14.095659672046486
preeq_numsteps      : [  0 463   0]
preeq_numsteps_b    : 0.0
preeq_status        : [<SteadyStateStatus.failed_factorization: -3>, <SteadyStateStatus.success: 1>, <SteadyStateStatus.not_run: 0>]
preeq_cpu_time      : 1.663
preeq_cpu_time_b    : 0.0

# Singular Jacobian, use simulation
model.set_steady_state_sensitivity_mode(
    SteadyStateSensitivityMode.integrateIfNewtonFails
)
solver = model.create_solver()
solver.set_newton_max_steps(10)
solver.set_sensitivity_method(SensitivityMethod.forward)
solver.set_sensitivity_order(SensitivityOrder.first)
rdata = run_simulation(model, solver, edata)
assert rdata.status == AMICI_SUCCESS
assert rdata.preeq_status == [
    SteadyStateStatus.failed_factorization,
    SteadyStateStatus.success,
    SteadyStateStatus.not_run,
]

for key, value in rdata.items():
    if key[0:6] == "preeq_":
        print(f"{key:20s}:", value)

preeq_wrms          : 1.1927149600038793e-08
preeq_t             : 32.0941543782181
preeq_numsteps      : [   0 1346    0]
preeq_numsteps_b    : 0.0
preeq_status        : [<SteadyStateStatus.failed_factorization: -3>, <SteadyStateStatus.success: 1>, <SteadyStateStatus.not_run: 0>]
preeq_cpu_time      : 6.43
preeq_cpu_time_b    : 0.0

# Non-singular Jacobian, use Newton solver
solver_reduced = model_reduced.create_solver()
solver_reduced.set_newton_max_steps(10)
solver_reduced.set_sensitivity_method(SensitivityMethod.forward)
solver_reduced.set_sensitivity_order(SensitivityOrder.first)
rdata_reduced = run_simulation(model_reduced, solver_reduced, edata)
assert rdata_reduced.status == AMICI_SUCCESS
for key, value in rdata_reduced.items():
    if key[0:6] == "preeq_":
        print(f"{key:20s}:", value)

preeq_wrms          : 7.777020191010815e-09
preeq_t             : nan
preeq_numsteps      : [2 0 0]
preeq_numsteps_b    : 0.0
preeq_status        : [<SteadyStateStatus.success: 1>, <SteadyStateStatus.not_run: 0>, <SteadyStateStatus.not_run: 0>]
preeq_cpu_time      : 0.024
preeq_cpu_time_b    : 0.0

Pre-equilibration with adjoint sensitivities

When using pre-equilibration, adjoint sensitivity analysis can be used for simulation. This is a particularly interesting case: Standard adjoint sensitivity analysis requires the initial state sensitivities sx0 to work, at least if data is given for finite (i.e., not exclusively post-equilibration) timepoints: For each parameter, a contribution to the gradient is given by the scalar product of the corresponding state sensitivity vector at timepoint \(t=t_0\), (column in sx0), with the adjoint state (\(p(t=t_0)\)). Hence, the matrix sx0 is needed. This scalar product “closes the loop” from forward to adjoint simulation.

By default, if adjoint sensitivity analysis is called with pre-equilibration, the initial state sensitivities are computed in just the same way as for forward sensitivity analysis. The only difference internally is that, if the steady state gets inferred via simulation, a separate solver object is used in order to ensure that the steady-state simulation does not interfere with the snapshotting of the forward trajectory from the actual time course.

However, also an adjoint version of pre-equilibration is possible: In this case, the “loop” from forward to adjoint simulation needs no closure: The simulation time is extended by pre-equilibration: forward from \(t = -\infty\) to \(t=t_0\), and after adjoint simulation also backward from \(t=t_0\) to \(t = -\infty\). Similar to adjoint post-equilibration, the steady state of the adjoint state (at \(t=-\infty\)) is \(p=0\), hence the scalar product (at \(t=-\infty\)) for the initial state sensitivities of pre-equilibration with the adjoint state vanishes. Instead, this gradient contribution is covered by additional quadratures \(\int_{-\infty}^0 p(s) ds \cdot \frac{\partial f}{\partial \theta}\). In order to compute these quadratures correctly, the adjoint state from the main adjoint simulation must be passed on to the initial adjoint state of backward pre-equilibration.

However, as the adjoint state must be passed on from backward computation to pre-equilibration, it is currently not allowed to alter (reinitialize) states of the model at \(t=t_0\), unless these states are constant, as otherwise this alteration would lead to a discontinuity in the adjoints state as well and hence to an incorrect gradient.

# Non-singular Jacobian, use Newton solver and adjoints with initial state sensitivities
solver_reduced = model_reduced.create_solver()
solver_reduced.set_newton_max_steps(10)
solver_reduced.set_sensitivity_method(SensitivityMethod.adjoint)
solver_reduced.set_sensitivity_order(SensitivityOrder.first)
rdata_reduced = run_simulation(model_reduced, solver_reduced, edata)

assert rdata_reduced.status == AMICI_SUCCESS
assert rdata_reduced.preeq_status == [
    SteadyStateStatus.success,
    SteadyStateStatus.not_run,
    SteadyStateStatus.not_run,
]
for key, value in rdata_reduced.items():
    if key[0:6] == "preeq_":
        print(f"{key:20s}:", value)
print("Gradient:", rdata_reduced["sllh"])

preeq_wrms          : 7.777020191010815e-09
preeq_t             : nan
preeq_numsteps      : [2 0 0]
preeq_numsteps_b    : 0.0
preeq_status        : [<SteadyStateStatus.success: 1>, <SteadyStateStatus.not_run: 0>, <SteadyStateStatus.not_run: 0>]
preeq_cpu_time      : 0.08
preeq_cpu_time_b    : 0.0
Gradient: [-2.33852861 -0.05841956  0.04908832 -0.03721072  6.33324895  0.        ]

# Non-singular Jacobian, use simulation solver and adjoints with initial state sensitivities
solver_reduced = model_reduced.create_solver()
solver_reduced.set_newton_max_steps(0)
solver_reduced.set_sensitivity_method(SensitivityMethod.adjoint)
solver_reduced.set_sensitivity_order(SensitivityOrder.first)
rdata_reduced = run_simulation(model_reduced, solver_reduced, edata)

assert rdata_reduced.status == AMICI_SUCCESS
assert rdata_reduced.preeq_status == [
    SteadyStateStatus.not_run,
    SteadyStateStatus.success,
    SteadyStateStatus.not_run,
]
assert rdata_reduced.preeq_numsteps_b == 0

for key, value in rdata_reduced.items():
    if key[0:6] == "preeq_":
        print(f"{key:20s}:", value)
print("Gradient:", rdata_reduced["sllh"])

preeq_wrms          : 0.9972941864795635
preeq_t             : 13.996121981443599
preeq_numsteps      : [  0 465   0]
preeq_numsteps_b    : 0.0
preeq_status        : [<SteadyStateStatus.not_run: 0>, <SteadyStateStatus.success: 1>, <SteadyStateStatus.not_run: 0>]
preeq_cpu_time      : 0.364
preeq_cpu_time_b    : 0.0
Gradient: [-2.33853656 -0.05841956  0.04908831 -0.03721076  6.33327232  0.        ]

# Non-singular Jacobian, use Newton solver and adjoints with fully adjoint pre-equilibration
solver_reduced = model_reduced.create_solver()
solver_reduced.set_newton_max_steps(10)
solver_reduced.set_sensitivity_method(SensitivityMethod.adjoint)
solver_reduced.set_sensitivity_method_pre_equilibration(
    SensitivityMethod.adjoint
)
solver_reduced.set_sensitivity_order(SensitivityOrder.first)
rdata_reduced = run_simulation(model_reduced, solver_reduced, edata)
assert rdata_reduced.status == AMICI_SUCCESS
assert rdata_reduced.preeq_status == [
    SteadyStateStatus.success,
    SteadyStateStatus.not_run,
    SteadyStateStatus.not_run,
]
assert rdata_reduced.preeq_numsteps_b == 0

for key, value in rdata_reduced.items():
    if key[0:6] == "preeq_":
        print(f"{key:20s}:", value)
print("Gradient:", rdata_reduced["sllh"])

preeq_wrms          : 7.777020191010815e-09
preeq_t             : nan
preeq_numsteps      : [2 0 0]
preeq_numsteps_b    : 0.0
preeq_status        : [<SteadyStateStatus.success: 1>, <SteadyStateStatus.not_run: 0>, <SteadyStateStatus.not_run: 0>]
preeq_cpu_time      : 0.017
preeq_cpu_time_b    : 0.006
Gradient: [-2.33909168 -0.05841956  0.04729539 -0.03889103  6.33324892  0.        ]

As for post-equilibration, adjoint pre-equilibration has an analytic solution (via the linear system), which will be preferred. If used for models with singular Jacobian, numerical integration will be carried out, which is indicated by preeq_numsteps_b.

# Singular Jacobian, use try Newton solver and adjoints with fully adjoint pre-equilibration
solver = model.create_solver()
solver.set_newton_max_steps(10)
solver.set_sensitivity_method(SensitivityMethod.adjoint)
solver.set_sensitivity_method_pre_equilibration(SensitivityMethod.adjoint)
solver.set_sensitivity_order(SensitivityOrder.first)
rdata = run_simulation(model, solver, edata)
assert rdata.status == AMICI_SUCCESS
assert rdata.preeq_status == [
    SteadyStateStatus.failed_factorization,
    SteadyStateStatus.success,
    SteadyStateStatus.not_run,
]

assert rdata.preeq_numsteps_b > 0
for key, value in rdata.items():
    if key[0:6] == "preeq_":
        print(f"{key:20s}:", value)
print("Gradient:", rdata["sllh"])

preeq_wrms          : 0.7841877852857544
preeq_t             : 14.095659672046486
preeq_numsteps      : [  0 463   0]
preeq_numsteps_b    : 605.0
preeq_status        : [<SteadyStateStatus.failed_factorization: -3>, <SteadyStateStatus.success: 1>, <SteadyStateStatus.not_run: 0>]
preeq_cpu_time      : 4.619
preeq_cpu_time_b    : 1.826
Gradient: [-2.32814211e+00 -5.84195575e-02  4.90883143e-02 -3.92894920e-02
  6.33327014e+00 -9.83788906e-09]

Controlling error tolerances in pre- and post-equilibration

When solving ODEs or DAEs, AMICI uses the default logic of CVODES and IDAS to control error tolerances. This means that error weights are computed based on the absolute error tolerances and the product of current state variables of the system and their respective relative error tolerances.

To find a steady state, numerical integration is performed until time derivatives \(\dot{\mathbf{x}}\) become sufficiently small. Namely, until the following condition is fulfilled:

\[\sqrt{\frac 1 n_x \sum_{i=1}^{n_x} (\dot{x}_i w_i)^2} < 1, \quad \text{where} \; w_i = \frac{1}{\text{rtol}*x_i + \text{atol}},\]

rtol and atol denote relative and absolute tolerances, respectively.

The respective tolerances for equilibrating a system with AMICI can be controlled by the user via the getter/setter functions [get|set][Absolute|Relative]ToleranceSteadyState[Sensi]:

# Non-singular Jacobian, use simulation
model_reduced.set_steady_state_sensitivity_mode(
    SteadyStateSensitivityMode.integrateIfNewtonFails
)
solver_reduced = model_reduced.create_solver()
solver_reduced.set_newton_max_steps(0)
solver_reduced.set_sensitivity_method(SensitivityMethod.forward)
solver_reduced.set_sensitivity_order(SensitivityOrder.first)

# run with lax tolerances
solver_reduced.set_relative_tolerance_steady_state(1e-2)
solver_reduced.set_absolute_tolerance_steady_state(1e-3)
solver_reduced.set_relative_tolerance_steady_state_sensi(1e-2)
solver_reduced.set_absolute_tolerance_steady_state_sensi(1e-3)
rdata_reduced_lax = run_simulation(model_reduced, solver_reduced, edata)

# run with strict tolerances
solver_reduced.set_relative_tolerance_steady_state(1e-12)
solver_reduced.set_absolute_tolerance_steady_state(1e-16)
solver_reduced.set_relative_tolerance_steady_state_sensi(1e-12)
solver_reduced.set_absolute_tolerance_steady_state_sensi(1e-16)
rdata_reduced_strict = run_simulation(model_reduced, solver_reduced, edata)

# compare ODE outputs
print("Number of ODE solver steps, until steady state:")
print("  lax tolerances:   ", rdata_reduced_lax["preeq_numsteps"])
print("  strict tolerances:", rdata_reduced_strict["preeq_numsteps"])

print("\nModel-internal time at which steady state was reached:")
print("  lax tolerances:   ", rdata_reduced_lax["preeq_t"])
print("  strict tolerances:", rdata_reduced_strict["preeq_t"])

print("\nCPU time to reach steady state (ms):")
print("  lax tolerances:   ", rdata_reduced_lax["preeq_cpu_time"])
print("  strict tolerances:", rdata_reduced_strict["preeq_cpu_time"])

Number of ODE solver steps, until steady state:
  lax tolerances:    [  0 881   0]
  strict tolerances: [   0 1262    0]

Model-internal time at which steady state was reached:
  lax tolerances:    6.741304044322333
  strict tolerances: 36.93904955664649

CPU time to reach steady state (ms):
  lax tolerances:    3.735
  strict tolerances: 5.282