Efficient computation of adjoint sensitivities at steady-state in ODE models of biochemical reaction networks

Dynamical models in the form of systems of ordinary differential equations have become a standard tool in systems biology. Many parameters of such models are usually unknown and have to be inferred from experimental data. Gradient-based optimization has proven to be effective for parameter estimation. However, computing gradients becomes increasingly costly for larger models, which are required for capturing the complex interactions of multiple biochemical pathways. Adjoint sensitivity analysis has been pivotal for working with such large models, but methods tailored for steady-state data are currently not available. We propose a new adjoint method for computing gradients, which is applicable if the experimental data include steady-state measurements. The method is based on a reformulation of the backward integration problem to a system of linear algebraic equations. The evaluation of the proposed method using real-world problems shows a speedup of total simulation time by a factor of up to 4.4. Our results demonstrate that the proposed approach can achieve a substantial improvement in computation time, in particular for large-scale models, where computational efficiency is critical.

We assume that the system 7 x = f (x(t, θ, u), θ, u), x(t 0 , θ) = x 0 (θ, u), has an exponentially stable steady state, which means that the real parts of the eigenvalues of the 8 Jacobian at the steady state (J(x * (θ, u), θ, u)) are negative. The system 9 p(t, θ, u) = −J(x * (θ, u), θ, u) T p(t, θ, u) has a steady state p = 0, which is then, as the real parts of the eigenvalues of J(x * (θ, u), θ, u) are negative, asymptotically stable in reverse time. Hence, on the interval [t ′ , t ′′ ), where the system (1) is at steady state, The computed p integral value is used to compute the objective function gradient by

10
In this section, we illustrate the proposed method on an exemplary conversion reaction where reaction rate coefficients θ 1 , θ 2 > 0. We perform the computations from Section "Methods,

11
Adjoint sensitivity analysis at steady state, Post-equilibration case" for this simple case.

S2
The ODE system describing these two reactions is The Jacobian of the system is which is singular; hence, the proposed method is not directly applicable.

13
In this case it is easy to see thatẋ A +ẋ B = 0 and the conserved quantity is x A (t) + x B (t). The 14 system can be simplified to where The Jacobian of this one-dimensional system is equal to −(θ 1 + θ 2 ) and does not depend on x A . This system has one non-trivial equilibrium which is exponentially stable as the only eigenvalue (−(θ 1 + θ 2 )) is negative.

16
When the system (3) is at steady state, the adjoint state is the solution oḟ and is equal to In this case The computed p integral value is used to compute the objective function gradient by S1.
The objective function gradient values were computed using ssASA for sensitivities and standard ASA for sensitivities. The values computed with the proposed ssASA method were considered accurate if the error between the two methods was small, with error (∆) calculated as where v ASA and v ssASA are the gradient values computed with the standard ASA or ssASA method for 19 sensitivities, respectively.

S3 AMICI implementation 21
The ASA at steady-state approach was implemented in the AMICI package. AMICI allows for forward 22 integration of differential equation models specified in SBML format or PySB, as well as for forward 23 sensitivity analysis, steady-state sensitivity analysis and ASA for likelihood-based output functions.

24
For ASA at steady-state, the user can choose between three options:

25
• only integration, which corresponds to numerical backward integration of the adjoint state ODE, 26 e.g. in the post-equilibration of the ODE on time intervals [t ′′ , t nt ), [t nt , t nt−1 ), . . . , [t 1 , t 0 ) with boundary values and lim t→t ′′+ p(t, θ, u) = 0 • only ssASA method, which corresponds to solving the linear system in the post-equilibration case, or system in pre-equilibration case.

28
• combined approach, where ssASA method is attempted first and, in case this fails, numerical 29 integration is used instead.

30
The alternative approaches are summarized in Fig S1.