{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Steady-state computation in AMICI\n", "\n", "*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.*\n", "\n", "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.\n", "\n", "AMICI offers several methods for computing steady states and sensitivities at steady state using both forward and adjoint sensitivity analysis.\n", "\n", "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.\n", "Common experimental scenarios include:\n", "\n", "1. the system is assumed to start in a steady state, then it is perturbed and enters a dynamic state;\n", "2. the system is assumed to start in a dynamic state, and after some time it reaches a steady state.\n", "\n", "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.\n", "\n", "![Simulation scenarios](figures/preeq_and_posteq.svg)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "> **_NOTE:_** Not every dynamical system needs to run into a steady state. Instead, it may exhibit:\n", ">\n", "> * continuous growth, e.g., $$\\dot{x} = x, \\quad x_0 = 1$$\n", "> * a finite-time blow up, e.g., $$\\dot{x} = x^2, \\quad x_0 = 1$$\n", "> * oscillations, e.g., $$\\ddot{x} = -x, \\quad x_0 = 1$$\n", "> * chaotic behaviour, e.g., the Lorenz attractor\n", "\n", "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.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The example model \n", "\n", "We will use the following model of an enzymatic reaction as example:\n", "\n", "![Model of an enzyme-catalyzed reaction](figures/catalyzed_reaction.svg)\n", "\n", "The corresponding system of ODEs is\n", "\n", "$$\n", "\\begin{array}{llll}\n", " \\dot{s} = (\\text{synthesis_substrate} + k_{\\text{create}}) - k_{\\text{bind}}\\cdot s \\cdot e + k_{\\text{unbind}}\\cdot c,\\\\\n", " \\dot{e} = - k_{\\text{bind}}\\cdot s \\cdot e + k_{\\text{unbind}}\\cdot c + k_{\\text{convert}}\\cdot c,\\\\\n", " \\dot{c} = k_{\\text{bind}}\\cdot s \\cdot e - k_{\\text{unbind}}\\cdot c - k_{\\text{convert}}\\cdot c,\\\\\n", " \\dot{p} = k_{\\text{convert}}\\cdot c - k_{\\text{decay}}\\cdot p,\n", "\\end{array}\n", "$$\n", "\n", "where $s$, $e$, $c$ and $p$ are concentrations of substrate, enzyme, complex and product, respectively.\n", "\n", "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\n", "\n", "$$\n", "\\begin{bmatrix}\n", " - k_{\\text{bind}}\\cdot e & - k_{\\text{bind}}\\cdot s & k_{\\text{unbind}} & 0\\\\\n", " - k_{\\text{bind}}\\cdot e & - k_{\\text{bind}}\\cdot s & (k_{\\text{unbind}} + k_{\\text{convert}}) & 0\\\\\n", " k_{\\text{bind}}\\cdot e & k_{\\text{bind}}\\cdot s & -(k_{\\text{unbind}} + k_{\\text{convert}}) & 0\\\\\n", " 0 & 0 & k_{\\text{convert}} & - k_{\\text{decay}}\\\\\n", "\\end{bmatrix}.\n", "$$\n", "\n", "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.\n", "\n", "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)](https://doi.org/10.1371/journal.pone.0100750).\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Import necessary libraries and define the model\n", "import os\n", "import tempfile\n", "from contextlib import suppress\n", "from pathlib import Path\n", "\n", "import amici\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from amici import (\n", " AMICI_ERROR,\n", " AMICI_SUCCESS,\n", " ExpData,\n", " SensitivityMethod,\n", " SensitivityOrder,\n", " SteadyStateComputationMode,\n", " SteadyStateSensitivityMode,\n", " SteadyStateStatus,\n", " import_model_module,\n", " run_simulation,\n", " simulation_status_to_str,\n", ")\n", "from amici import (\n", " MeasurementChannel as MC,\n", ")\n", "from amici.antimony_import import antimony2sbml\n", "from amici.plotting import (\n", " plot_observable_trajectories,\n", " plot_state_trajectories,\n", ")\n", "\n", "# We encode the model in Antimony format (https://tellurium.readthedocs.io/en/latest/antimony.html),\n", "# which is a human-readable format for biochemical models that can be converted to SBML.\n", "ant_model = \"\"\"model *model_constant_species()\n", " const compartment compartment_ = 1;\n", "\n", " // species\n", " var species substrate in compartment_ = 0;\n", " var species enzyme in compartment_;\n", " var species complex in compartment_ = 0;\n", " species product in compartment_ = 0;\n", "\n", " // dynamic expressions\n", " var substrate_product := complex + product + substrate;\n", "\n", " // reactions\n", " creation: => substrate; compartment_*(synthesis_substrate + k_create);\n", " binding: substrate + enzyme -> complex; compartment_*(k_bind*substrate*enzyme - k_unbind*complex);\n", " conversion: complex => enzyme + product; compartment_*k_convert*complex;\n", " decay: product => ; compartment_*k_decay*product;\n", "\n", " // initial assignments\n", " enzyme = init_enzyme;\n", " substrate = init_substrate;\n", "\n", " // parameters\n", " const init_enzyme = 2;\n", " const init_substrate = 5;\n", " const synthesis_substrate = 0;\n", " const k_bind = 10;\n", " const k_unbind = 1;\n", " const k_convert = 10;\n", " const k_create = 2;\n", " const k_decay = 1;\n", "\n", " unit volume = 1e-3 litre;\n", " unit substance = 1e-3 mole;\n", "end\"\"\"\n", "\n", "\n", "# Name of the models that will also be the name of the python module\n", "model_name = \"model_constant_species\"\n", "model_reduced_name = model_name + \"_reduced\"\n", "\n", "# Directories to which the generated model code is written\n", "temp_dir = Path(tempfile.mkdtemp())\n", "model_output_dir = temp_dir / model_name\n", "model_reduced_output_dir = temp_dir / model_reduced_name" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Import the model\n", "sbml_importer = amici.SbmlImporter(antimony2sbml(ant_model), from_file=False)\n", "\n", "# specify observables and constant parameters\n", "constant_parameters = [\"synthesis_substrate\", \"init_enzyme\"]\n", "observation_model = [\n", " MC(id_=\"observable_product\", name=\"\", formula=\"product\"),\n", " MC(id_=\"observable_substrate\", name=\"\", formula=\"substrate\"),\n", "]\n", "sigmas = {\"observable_product\": 1.0, \"observable_substrate\": 1.0}\n", "\n", "# import the model without handled the conserved quantity\n", "with suppress(KeyError):\n", " del os.environ[\"AMICI_EXPERIMENTAL_SBML_NONCONST_CLS\"]\n", "sbml_importer.sbml2amici(\n", " model_name,\n", " model_output_dir,\n", " observation_model=observation_model,\n", " constant_parameters=constant_parameters,\n", " compute_conservation_laws=False,\n", ")\n", "\n", "# turn on removal of conservation laws and import the reduced model\n", "os.environ[\"AMICI_EXPERIMENTAL_SBML_NONCONST_CLS\"] = \"1\"\n", "sbml_importer.sbml2amici(\n", " model_reduced_name,\n", " model_reduced_output_dir,\n", " observation_model=observation_model,\n", " constant_parameters=constant_parameters,\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# import the models and run some test simulations\n", "model_reduced_module = import_model_module(\n", " model_reduced_name, model_reduced_output_dir\n", ")\n", "model_reduced = model_reduced_module.get_model()\n", "\n", "model_module = import_model_module(model_name, model_output_dir)\n", "model = model_module.get_model()\n", "\n", "# simulate model with conservation laws\n", "model.set_timepoints(np.linspace(0, 2, 100))\n", "solver = model.create_solver()\n", "rdata = run_simulation(model, solver)\n", "assert rdata.status == AMICI_SUCCESS\n", "\n", "# simulate model without conservation laws\n", "model_reduced.set_timepoints(np.linspace(0, 2, 100))\n", "solver_reduced = model_reduced.create_solver()\n", "rdata_reduced = run_simulation(model_reduced, solver_reduced)\n", "assert rdata_reduced.status == AMICI_SUCCESS\n", "\n", "# plot trajectories\n", "fig, axes = plt.subplots(2, 2, figsize=(8, 8))\n", "plot_state_trajectories(rdata, model=model, ax=axes[0, 0])\n", "plot_observable_trajectories(rdata, model=model, ax=axes[1, 0])\n", "plot_state_trajectories(rdata_reduced, model=model_reduced, ax=axes[0, 1])\n", "plot_observable_trajectories(rdata_reduced, model=model_reduced, ax=axes[1, 1])\n", "fig.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# the enzyme state was removed from the ODEs of the reduced model\n", "print(model.get_state_ids_solver())\n", "print(model_reduced.get_state_ids_solver())\n", "print(model.get_state_ids())\n", "print(model_reduced.get_state_ids())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Steady-state computation \n", "\n", "Two approaches are available in AMICI for steady-state computation:\n", "\n", "1. Newton's method\n", "\n", "2. Numerical integration until time derivatives $\\dot{x}$ become sufficiently small.\n", " Namely, numerical integration is performed until the condition\n", "\n", " $$\n", " \\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}}\n", " $$\n", "\n", " is fulfilled, where “rtol” and “atol” denote relative and absolute tolerances, respectively.\n", "\n", "Newton's method can be significantly faster than numerical integration.\n", "However, it will converge only if it was started close enough to the steady state.\n", "More about the convergence of the Newton's method and its implementation in AMICI can be found in [Lines et al. (2019)](https://doi.org/10.1016/j.ifacol.2019.12.232).\n", "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)](https://doi.org/10.1371/journal.pone.0312148)).\n", " Therefore, AMICI users numerical integration by default, which is more robust and can be used for any model.\n", "\n", "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.\n", "\n", "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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Inferring the steady state of the system (post-equilibration) \n", "\n", "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.\n", "First, we want to demonstrate that Newton's method will fail with the unreduced model due to a singular right-hand side Jacobian." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Call post-equilibration by setting an infinity timepoint\n", "model.set_timepoints([np.inf])\n", "\n", "# Set the solver\n", "solver = model.create_solver()\n", "solver.set_newton_max_steps(10)\n", "# maximum number of solver steps for numerical integration\n", "solver.set_max_steps(1000)\n", "model.set_steady_state_computation_mode(\n", " SteadyStateComputationMode.integrateIfNewtonFails\n", ")\n", "rdata = run_simulation(model, solver)\n", "assert rdata.status == AMICI_SUCCESS\n", "# print out a subset of data returned by model simulation\n", "for key in [\n", " \"ts\",\n", " \"x\",\n", " \"xdot\",\n", " \"posteq_wrms\",\n", " \"posteq_t\",\n", " \"posteq_numsteps\",\n", " \"posteq_numsteps_b\",\n", " \"posteq_status\",\n", " \"posteq_cpu_time\",\n", " \"posteq_cpu_time_b\",\n", "]:\n", " print(f\"{key:>16}: {rdata[key]}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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`.\n", "\n", "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\n", "\n", " 1. the Newton's method run;\n", " 2. the numerical integration run;\n", " 3. the second Newton's method run, starting from the simulation results.\n", "\n", "The status is of type `SteadyStateStatus` with the following meanings:\n", "\n", " * `SteadyStateStatus.success`: Successful run\n", " * `SteadyStateStatus.not_run`: Did not run\n", " * `SteadyStateStatus.failed`: Error: No further specification is given, the error message should give more information.\n", " * `SteadyStateStatus.failed_convergence`: Error: The method did not converge to a steady state within the maximum number of steps (Newton's method or simulation).\n", " * `SteadyStateStatus.failed_factorization`: Error: The Jacobian of the right hand side is singular (only Newton's method)\n", " * `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)\n", " * `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.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "list(SteadyStateStatus)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In our case, only the second entry of `posteq_status` contains a positive integer:\n", "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).\n", "The third entry is `0`, thus Newton's method was not launched for the second time.\n", "More information can be found in`posteq_numsteps`:\n", "Also here, only the second entry contains a positive integer, which is smaller than the allowed maximum number of steps (1000).\n", "Hence, steady state was reached via simulation, which corresponds to the simulated time written to `posteq_time`.\n", "\n", "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]`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# steady state value found by pre-equilibration\n", "steady_state = rdata[\"x\"][0]\n", "\n", "timepoints = np.linspace(0, 14, 200)\n", "model.set_timepoints(timepoints)\n", "rdata = run_simulation(model, solver)\n", "plot_state_trajectories(rdata, model=model)\n", "\n", "for stst_value in steady_state:\n", " plt.axhline(y=stst_value, color=\"gray\", linestyle=\"--\", linewidth=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": "We want to demonstrate a complete failure during the steady-state computation by reducing the number of integration steps to a lower value:" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# reduce maxsteps for integration\n", "model.set_timepoints([np.inf])\n", "solver.set_max_steps(100)\n", "rdata = run_simulation(model, solver)\n", "assert rdata.status == AMICI_ERROR\n", "print(\"Simulation status:\", simulation_status_to_str(rdata[\"status\"]))\n", "print(\"Status of post-equilibration:\", rdata[\"posteq_status\"])\n", "print(\n", " \"Number of steps employed in post-equilibration:\", rdata[\"posteq_numsteps\"]\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, the same logic works if we use the reduced model.\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model_reduced.set_timepoints([np.inf])\n", "model_reduced.set_steady_state_computation_mode(\n", " SteadyStateComputationMode.integrateIfNewtonFails\n", ")\n", "\n", "# set the solver\n", "solver_reduced = model_reduced.create_solver()\n", "solver_reduced.set_newton_max_steps(10)\n", "solver_reduced.set_max_steps(100)\n", "rdata_reduced = run_simulation(model_reduced, solver_reduced)\n", "assert rdata_reduced.status == AMICI_SUCCESS\n", "assert rdata_reduced.posteq_status[0] == SteadyStateStatus.success\n", "print(\n", " \"Simulation status:\",\n", " simulation_status_to_str(rdata_reduced[\"status\"]),\n", ")\n", "print(\"Steady-state value:\", rdata_reduced[\"x\"][0])\n", "print(\"Status of post-equilibration:\", rdata_reduced[\"posteq_status\"])\n", "print(\n", " \"Number of steps employed in post-equilibration:\",\n", " rdata_reduced[\"posteq_numsteps\"],\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Pre-equilibrating the model\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# create edata, with 3 timepoints and 2 observables:\n", "edata = ExpData(2, 0, 0, np.array([0.0, 0.1, 1.0]))\n", "edata.set_observed_data([1.8] * 6)\n", "edata.fixed_parameters = np.array([3.0, 5.0])\n", "# set parameters for pre-equilibration\n", "edata.fixed_parameters_pre_equilibration = np.array([0.0, 2.0])\n", "edata.reinitialize_fixed_parameter_initial_states = True" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# create the solver object and run the simulation\n", "solver_reduced = model_reduced.create_solver()\n", "solver_reduced.set_newton_max_steps(10)\n", "rdata_reduced = run_simulation(model_reduced, solver_reduced, edata)\n", "assert rdata_reduced.status == AMICI_SUCCESS\n", "print(\n", " \"Simulation status:\",\n", " simulation_status_to_str(rdata_reduced[\"status\"]),\n", ")\n", "print(\"Status of pre-equilibration:\", rdata_reduced[\"preeq_status\"])\n", "print(\n", " \"Number of steps employed in pre-equilibration:\",\n", " rdata_reduced[\"preeq_numsteps\"],\n", ")\n", "print(\"Steady state found during pre-equilibration:\", rdata_reduced[\"x_ss\"])\n", "print(\"Initial state used for simulation:\", rdata_reduced[\"x0\"])\n", "\n", "\n", "plot_state_trajectories(rdata_reduced, model=model_reduced)\n", "plot_observable_trajectories(rdata_reduced, model=model_reduced)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "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.\n", "\n", "We can also combine pre- and post-equilibration:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Change the last timepoint to an infinity timepoint.\n", "edata.set_timepoints(np.array([0.0, 0.1, float(\"inf\")]))\n", "\n", "# run the simulation\n", "rdata_reduced = run_simulation(model_reduced, solver_reduced, edata)\n", "assert rdata_reduced.status == AMICI_SUCCESS" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sensitivities computation at steady state \n", "\n", "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.\n", "\n", "### Post-equilibration with sensitivities \n", "\n", "\n", "\n", "#### Post-equilibration with forward sensitivities\n", "\n", "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:\n", "\n", "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.\n", "\n", "\n", "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:\n", " $$\n", " 0 = \\left.\\dot{\\mathbf{s}}^x\\right\\vert_{\\mathbf{x} = \\mathbf{x}^*({\\theta})}\n", " = \\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})}\n", " \\qquad \\Rightarrow \\qquad\n", " \\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})},$$\n", " where $\\mathbf{s}^x$ are the state sensitivities,\n", " $\\mathbf{f}$ is the right-hand side of the model ODEs,\n", " $\\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.\n", " However, this method is only applicable if the Jacobian is not (close to) singular.\n", "\n", "\n", "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).\n", "\n", "\n", "Side remark:\n", "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.\n", "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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "\n", "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Call simulation with singular Jacobian and `integrateIfNewtonFails` mode\n", "model.set_timepoints([np.inf])\n", "model.set_steady_state_sensitivity_mode(\n", " SteadyStateSensitivityMode.integrateIfNewtonFails\n", ")\n", "solver = model.create_solver()\n", "solver.set_newton_max_steps(10)\n", "solver.set_sensitivity_method(SensitivityMethod.forward)\n", "solver.set_sensitivity_order(SensitivityOrder.first)\n", "solver.set_max_steps(10000)\n", "rdata = run_simulation(model, solver)\n", "\n", "np.set_printoptions(threshold=20)\n", "print(\"Simulation status:\", simulation_status_to_str(rdata[\"status\"]))\n", "print(\"Status of post-equilibration:\", rdata[\"posteq_status\"])\n", "print(\n", " \"Number of steps employed in post-equilibration:\", rdata[\"posteq_numsteps\"]\n", ")\n", "print(\"Computed state sensitivities:\")\n", "print(rdata[\"sx\"][0, :, :])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "The shape of `rdata['sx']` is (number of timepoints) x (number of parameters) x (number of states)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Call simulation with singular Jacobian and newtonOnly mode (will fail)\n", "model.set_timepoints([np.inf])\n", "model.set_steady_state_sensitivity_mode(SteadyStateSensitivityMode.newtonOnly)\n", "solver = model.create_solver()\n", "solver.set_newton_max_steps(10)\n", "solver.set_sensitivity_method(SensitivityMethod.forward)\n", "solver.set_sensitivity_order(SensitivityOrder.first)\n", "solver.set_max_steps(10000)\n", "rdata = run_simulation(model, solver)\n", "assert rdata.status == AMICI_ERROR\n", "\n", "print(\"Simulation status:\", simulation_status_to_str(rdata[\"status\"]))\n", "print(\"Steady state:\", rdata[\"x\"])\n", "print(\"Status of post-equilibration:\", rdata[\"posteq_status\"])\n", "print(\n", " \"Number of steps employed in post-equilibration:\", rdata[\"posteq_numsteps\"]\n", ")\n", "print(\"Computed state sensitivities:\")\n", "print(rdata[\"sx\"][0, :, :])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the model with singular Jacobian steady-state computation worked by numerical integration.\n", "However, for sensitivities computation `SteadyStateSensitivityMode.newtonOnly` approach was chosen which fails with singular Jacobian and the sensitivities could not be computed.\n", "We can try the same with the reduced model:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Try `newtonOnly` option with reduced model\n", "model_reduced.set_timepoints([np.inf])\n", "model_reduced.set_steady_state_sensitivity_mode(\n", " SteadyStateSensitivityMode.newtonOnly\n", ")\n", "solver_reduced = model_reduced.create_solver()\n", "solver_reduced.set_newton_max_steps(10)\n", "solver_reduced.set_sensitivity_method(SensitivityMethod.forward)\n", "solver_reduced.set_sensitivity_order(SensitivityOrder.first)\n", "solver_reduced.set_max_steps(1000)\n", "rdata_reduced = run_simulation(model_reduced, solver_reduced)\n", "assert rdata_reduced.status == AMICI_SUCCESS\n", "assert rdata_reduced.posteq_status == [\n", " SteadyStateStatus.success,\n", " SteadyStateStatus.not_run,\n", " SteadyStateStatus.not_run,\n", "]\n", "\n", "print(\n", " \"Simulation status:\",\n", " simulation_status_to_str(rdata_reduced[\"status\"]),\n", ")\n", "print(\"Steady state:\", rdata_reduced[\"x\"])\n", "print(\"Status of post-equilibration:\", rdata_reduced[\"posteq_status\"])\n", "print(\n", " \"Number of steps employed in post-equilibration:\",\n", " rdata_reduced[\"posteq_numsteps\"],\n", ")\n", "print(\"Computed state sensitivities:\")\n", "print(rdata_reduced[\"sx\"][0, :, :])" ] }, { "cell_type": "markdown", "metadata": {}, "source": "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$)." }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Post-equilibration with adjoint sensitivities\n", "\n", "Post-equilibration also works with adjoint sensitivity analysis and, similar to forward sensitivities, two approaches are available for sensitivities computation at steady state:\n", "\n", "1. Standard adjoint sensitivity analysis can be performed to compute the objective function gradient.\n", " 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}$).\n", " Ultimately, one needs to compute\n", "$$\n", "\\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,\n", "$$\n", "where $\\mathbf{p}$ is the adjoint state, $\\theta$ the model parameters.\n", "\n", "\n", "2. Alternatively, one can exploit the steady-state condition and speed up the computation.\n", " It can be shown that the integral\n", "$$\n", "\\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\n", "$$\n", " can be reduced to a matrix-vector product\n", "$$\n", "\\int^{t''}_{t'}\n", "\\mathbf{p}(t, {\\theta}, \\mathbf{u})^T dt \\cdot \\left.\n", "\\frac{\\partial \\mathbf{f}}\n", "{\\partial \\theta_k}\n", "\\right\\vert_{\\mathbf{x} = \\mathbf{x}^*({\\theta})}dt\n", "=\n", "\\mathbf{p}_{\\text{integral}}\n", "\\cdot \\left.\n", "\\frac{\\partial \\mathbf{f}}\n", "{\\partial \\theta_k}\n", "\\right\\vert_{\\mathbf{x} = \\mathbf{x}^*({\\theta})},\n", "$$\n", " 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\n", "$$\n", " \\mathbf{J}(\\mathbf{x}^*({\\theta}), {\\theta})^T \\mathbf{p}_{\\text{integral}} = -\\mathbf{p}(t'', {\\theta}).\n", "$$\n", " The detailed explanation of this approach and the derivation of the last equation can be found in [Lakrisenko et al. (2023)](https://doi.org/10.1371/journal.pcbi.1010783).\n", "\n", " 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.\n", "\n", "![Adjoint steady-state sensitivity computation modes in AMICI](figures/amici_logic_adjoint_sensitivities_at_stst.svg)\n", "\n", "**_NOTE:_** adjoint sensitivity analysis yields less information than forward sensitivity analysis,\n", "since state sensitivities are not computed. However, it has a better scaling behavior towards large model size.\n", "\n", "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Call post-equilibration and sensitivities computation using adjoint sensitivity analysis\n", "# by setting an infinity timepoint\n", "# and creatíng an edata object, which is needed for adjoint computation\n", "edata = ExpData(2, 0, 0, np.array([float(\"inf\")]))\n", "edata.set_observed_data([1.8] * 2)\n", "edata.fixed_parameters = np.array([3.0, 5.0])\n", "\n", "model_reduced.set_steady_state_sensitivity_mode(\n", " SteadyStateSensitivityMode.newtonOnly\n", ")\n", "solver_reduced = model_reduced.create_solver()\n", "solver_reduced.set_newton_max_steps(10)\n", "solver_reduced.set_sensitivity_method(SensitivityMethod.adjoint)\n", "solver_reduced.set_sensitivity_order(SensitivityOrder.first)\n", "solver_reduced.set_max_steps(1000)\n", "rdata_reduced = run_simulation(model_reduced, solver_reduced, edata)\n", "\n", "print(\n", " \"Simulation status:\",\n", " simulation_status_to_str(rdata_reduced[\"status\"]),\n", ")\n", "print(\"Status of post-equilibration:\", rdata_reduced[\"posteq_status\"])\n", "print(\n", " \"Number of steps employed in post-equilibration:\",\n", " rdata_reduced[\"posteq_numsteps\"],\n", ")\n", "print(\n", " \"Number of backward steps employed in post-equilibration:\",\n", " rdata_reduced[\"posteq_numsteps_b\"],\n", ")\n", "print(\"Computed gradient:\", rdata_reduced[\"sllh\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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).\n", "Now, integration is carried out and hence `posteq_numsteps_b > 0`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Call adjoint post-equilibration with model with singular Jacobian\n", "model.set_steady_state_sensitivity_mode(SteadyStateSensitivityMode.newtonOnly)\n", "solver = model.create_solver()\n", "solver.set_newton_max_steps(10)\n", "solver.set_max_steps_backward_problem(10000)\n", "solver.set_sensitivity_method(SensitivityMethod.adjoint)\n", "solver.set_sensitivity_order(SensitivityOrder.first)\n", "rdata = run_simulation(model, solver, edata)\n", "\n", "print(\"Simulation status:\", simulation_status_to_str(rdata[\"status\"]))\n", "print(\"Status of post-equilibration:\", rdata[\"posteq_status\"])\n", "print(\n", " \"Number of steps employed in post-equilibration:\", rdata[\"posteq_numsteps\"]\n", ")\n", "print(\n", " \"Number of backward steps employed in post-equilibration:\",\n", " rdata[\"posteq_numsteps_b\"],\n", ")\n", "print(\"Computed gradient:\", rdata[\"sllh\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Pre-equilibration with sensitivities \n", "\n", "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.\n", "However, if forward simulation with steady state sensitivities is allowed, or if the Jacobian is not singular, it will work." ] }, { "cell_type": "markdown", "metadata": {}, "source": "#### Pre-equilibration with forward sensitivities" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# No post-equilibration this time.\n", "# create edata, with 3 timepoints and 2 observables:\n", "edata = ExpData(2, 0, 0, np.array([0.0, 0.1, 1.0]))\n", "edata.set_observed_data([1.8] * 6)\n", "edata.fixed_parameters = np.array([3.0, 5.0])\n", "# set parameters for pre-equilibration\n", "edata.fixed_parameters_pre_equilibration = np.array([0.0, 2.0])\n", "edata.reinitialize_fixed_parameter_initial_states = True\n", "\n", "# create the solver object and run the simulation, singular Jacobian, enforce Newton solver for sensitivities\n", "model.set_steady_state_sensitivity_mode(SteadyStateSensitivityMode.newtonOnly)\n", "solver = model.create_solver()\n", "solver.set_newton_max_steps(10)\n", "solver.set_sensitivity_method(SensitivityMethod.forward)\n", "solver.set_sensitivity_order(SensitivityOrder.first)\n", "rdata = run_simulation(model, solver, edata)\n", "\n", "assert rdata.status == AMICI_ERROR\n", "for key, value in rdata.items():\n", " if key[0:6] == \"preeq_\":\n", " print(f\"{key:20s}:\", value)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Singular Jacobian, use simulation\n", "model.set_steady_state_sensitivity_mode(\n", " SteadyStateSensitivityMode.integrateIfNewtonFails\n", ")\n", "solver = model.create_solver()\n", "solver.set_newton_max_steps(10)\n", "solver.set_sensitivity_method(SensitivityMethod.forward)\n", "solver.set_sensitivity_order(SensitivityOrder.first)\n", "rdata = run_simulation(model, solver, edata)\n", "assert rdata.status == AMICI_SUCCESS\n", "assert rdata.preeq_status == [\n", " SteadyStateStatus.failed_factorization,\n", " SteadyStateStatus.success,\n", " SteadyStateStatus.not_run,\n", "]\n", "\n", "for key, value in rdata.items():\n", " if key[0:6] == \"preeq_\":\n", " print(f\"{key:20s}:\", value)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Non-singular Jacobian, use Newton solver\n", "solver_reduced = model_reduced.create_solver()\n", "solver_reduced.set_newton_max_steps(10)\n", "solver_reduced.set_sensitivity_method(SensitivityMethod.forward)\n", "solver_reduced.set_sensitivity_order(SensitivityOrder.first)\n", "rdata_reduced = run_simulation(model_reduced, solver_reduced, edata)\n", "assert rdata_reduced.status == AMICI_SUCCESS\n", "for key, value in rdata_reduced.items():\n", " if key[0:6] == \"preeq_\":\n", " print(f\"{key:20s}:\", value)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Pre-equilibration with adjoint sensitivities\n", "\n", "\n", "When using pre-equilibration, adjoint sensitivity analysis can be used for simulation.\n", "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:\n", "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.\n", "\n", "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.\n", "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.\n", "\n", "However, also an adjoint version of pre-equilibration is possible:\n", "In this case, the \"loop\" from forward to adjoint simulation needs no closure:\n", "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$.\n", "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.\n", "Instead, this gradient contribution is covered by additional quadratures $\\int_{-\\infty}^0 p(s) ds \\cdot \\frac{\\partial f}{\\partial \\theta}$.\n", "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.\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Non-singular Jacobian, use Newton solver and adjoints with initial state sensitivities\n", "solver_reduced = model_reduced.create_solver()\n", "solver_reduced.set_newton_max_steps(10)\n", "solver_reduced.set_sensitivity_method(SensitivityMethod.adjoint)\n", "solver_reduced.set_sensitivity_order(SensitivityOrder.first)\n", "rdata_reduced = run_simulation(model_reduced, solver_reduced, edata)\n", "\n", "assert rdata_reduced.status == AMICI_SUCCESS\n", "assert rdata_reduced.preeq_status == [\n", " SteadyStateStatus.success,\n", " SteadyStateStatus.not_run,\n", " SteadyStateStatus.not_run,\n", "]\n", "for key, value in rdata_reduced.items():\n", " if key[0:6] == \"preeq_\":\n", " print(f\"{key:20s}:\", value)\n", "print(\"Gradient:\", rdata_reduced[\"sllh\"])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Non-singular Jacobian, use simulation solver and adjoints with initial state sensitivities\n", "solver_reduced = model_reduced.create_solver()\n", "solver_reduced.set_newton_max_steps(0)\n", "solver_reduced.set_sensitivity_method(SensitivityMethod.adjoint)\n", "solver_reduced.set_sensitivity_order(SensitivityOrder.first)\n", "rdata_reduced = run_simulation(model_reduced, solver_reduced, edata)\n", "\n", "assert rdata_reduced.status == AMICI_SUCCESS\n", "assert rdata_reduced.preeq_status == [\n", " SteadyStateStatus.not_run,\n", " SteadyStateStatus.success,\n", " SteadyStateStatus.not_run,\n", "]\n", "assert rdata_reduced.preeq_numsteps_b == 0\n", "\n", "for key, value in rdata_reduced.items():\n", " if key[0:6] == \"preeq_\":\n", " print(f\"{key:20s}:\", value)\n", "print(\"Gradient:\", rdata_reduced[\"sllh\"])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Non-singular Jacobian, use Newton solver and adjoints with fully adjoint pre-equilibration\n", "solver_reduced = model_reduced.create_solver()\n", "solver_reduced.set_newton_max_steps(10)\n", "solver_reduced.set_sensitivity_method(SensitivityMethod.adjoint)\n", "solver_reduced.set_sensitivity_method_pre_equilibration(\n", " SensitivityMethod.adjoint\n", ")\n", "solver_reduced.set_sensitivity_order(SensitivityOrder.first)\n", "rdata_reduced = run_simulation(model_reduced, solver_reduced, edata)\n", "assert rdata_reduced.status == AMICI_SUCCESS\n", "assert rdata_reduced.preeq_status == [\n", " SteadyStateStatus.success,\n", " SteadyStateStatus.not_run,\n", " SteadyStateStatus.not_run,\n", "]\n", "assert rdata_reduced.preeq_numsteps_b == 0\n", "\n", "for key, value in rdata_reduced.items():\n", " if key[0:6] == \"preeq_\":\n", " print(f\"{key:20s}:\", value)\n", "print(\"Gradient:\", rdata_reduced[\"sllh\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": "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`." }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Singular Jacobian, use try Newton solver and adjoints with fully adjoint pre-equilibration\n", "solver = model.create_solver()\n", "solver.set_newton_max_steps(10)\n", "solver.set_sensitivity_method(SensitivityMethod.adjoint)\n", "solver.set_sensitivity_method_pre_equilibration(SensitivityMethod.adjoint)\n", "solver.set_sensitivity_order(SensitivityOrder.first)\n", "rdata = run_simulation(model, solver, edata)\n", "assert rdata.status == AMICI_SUCCESS\n", "assert rdata.preeq_status == [\n", " SteadyStateStatus.failed_factorization,\n", " SteadyStateStatus.success,\n", " SteadyStateStatus.not_run,\n", "]\n", "\n", "assert rdata.preeq_numsteps_b > 0\n", "for key, value in rdata.items():\n", " if key[0:6] == \"preeq_\":\n", " print(f\"{key:20s}:\", value)\n", "print(\"Gradient:\", rdata[\"sllh\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Controlling error tolerances in pre- and post-equilibration\n", "\n", "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.\n", "\n", "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:\n", "\n", "$$\n", " \\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}},\n", "$$\n", "`rtol` and `atol` denote relative and absolute tolerances, respectively.\n", "\n", "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]`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Non-singular Jacobian, use simulation\n", "model_reduced.set_steady_state_sensitivity_mode(\n", " SteadyStateSensitivityMode.integrateIfNewtonFails\n", ")\n", "solver_reduced = model_reduced.create_solver()\n", "solver_reduced.set_newton_max_steps(0)\n", "solver_reduced.set_sensitivity_method(SensitivityMethod.forward)\n", "solver_reduced.set_sensitivity_order(SensitivityOrder.first)\n", "\n", "# run with lax tolerances\n", "solver_reduced.set_relative_tolerance_steady_state(1e-2)\n", "solver_reduced.set_absolute_tolerance_steady_state(1e-3)\n", "solver_reduced.set_relative_tolerance_steady_state_sensi(1e-2)\n", "solver_reduced.set_absolute_tolerance_steady_state_sensi(1e-3)\n", "rdata_reduced_lax = run_simulation(model_reduced, solver_reduced, edata)\n", "\n", "# run with strict tolerances\n", "solver_reduced.set_relative_tolerance_steady_state(1e-12)\n", "solver_reduced.set_absolute_tolerance_steady_state(1e-16)\n", "solver_reduced.set_relative_tolerance_steady_state_sensi(1e-12)\n", "solver_reduced.set_absolute_tolerance_steady_state_sensi(1e-16)\n", "rdata_reduced_strict = run_simulation(model_reduced, solver_reduced, edata)\n", "\n", "# compare ODE outputs\n", "print(\"Number of ODE solver steps, until steady state:\")\n", "print(\" lax tolerances: \", rdata_reduced_lax[\"preeq_numsteps\"])\n", "print(\" strict tolerances:\", rdata_reduced_strict[\"preeq_numsteps\"])\n", "\n", "print(\"\\nModel-internal time at which steady state was reached:\")\n", "print(\" lax tolerances: \", rdata_reduced_lax[\"preeq_t\"])\n", "print(\" strict tolerances:\", rdata_reduced_strict[\"preeq_t\"])\n", "\n", "print(\"\\nCPU time to reach steady state (ms):\")\n", "print(\" lax tolerances: \", rdata_reduced_lax[\"preeq_cpu_time\"])\n", "print(\" strict tolerances:\", rdata_reduced_strict[\"preeq_cpu_time\"])" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.6" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }