Balanced truncation for model reduction of biological oscillators

Model reduction is a central problem in mathematical biology. Reduced order models enable modeling of a biological system at different levels of complexity and the quantitative analysis of its properties, like sensitivity to parameter variations and resilience to exogenous perturbations. However, available model reduction methods often fail to capture a diverse range of nonlinear behaviors observed in biology, such as multistability and limit cycle oscillations. The paper addresses this need using differential analysis. This approach leads to a nonlinear enhancement of classical balanced truncation for biological systems whose behavior is not restricted to the stability of a single equilibrium. Numerical results suggest that the proposed framework may be relevant to the approximation of classical models of biological systems.


Introduction
Model reduction is a central problem in mathematical biology (Fall et al. 2005;Murray 2007;Alon 2007;Keener and Sneyd 2008;Del Vecchio and Murray 2015). Besides their role in numerical simulations, models enable the study of the principles underlying the behavior of a biological system and the quantitative analysis of its properties, like sensitivity to parameter variations and resilience to exogenous perturbations. Models of biological phenomena often originate from complex networks of chemical reactions (Fall et al. 2005;Murray 2007;Alon 2007;Keener and Sneyd 2008;Del Vecchio and Murray 2015), the temporal dynamics of which is modeled by possibly large systems of differential equations. The resulting nonlinear dynamical models pose significant Communicated by Jordi Garcia-Ojalvo. challenges for simulation, analysis and design. Model reduction alleviates these issues by constructing reduced order models whose behavior captures that of the original system (Antoulas 2005). Motivated by the goal of designing biological devices of increasing size and complexity , model reduction of biological systems has recently attracted renewed interest in the rapidly emerging field of systems and synthetic biology (Gómez-Uribe et al. 2008;Anderson et al. 2011;Thomas et al. 2012;Kang and Kurtz 2013;Rao et al. 2014;Sootla and Anderson 2014;Radulescu et al. 2015;Del Vecchio and Murray 2015;Herath et al. 2016;Sootla and Anderson 2017;Herath and Del Vecchio 2018). Given the pressing need for compositional modeling frameworks in mathematical biology , the present paper seeks to develop a model reduction framework compatible with modeling and interconnection of open systems whose behavior is not restricted to the stability of a single equilibrium.
Background Model reduction of biological systems has a long history (Snowden et al. 2017), with early contributions dating back as far as the 1920s (Briggs and Haldane 1925). The classical approach to model reduction of biological systems relies on timescale separation arguments (Briggs and Haldane 1925;Segel and Slemrod 1989;Gómez-Uribe et al. 2008; Thomas et al. 2012;Kang and Kurtz 2013;Prescott and Papachristodoulou 2014;Del Vecchio and Murray 2015;Herath et al. 2016;Herath and Del Vecchio 2018). In general, many biochemical processes are described by reactions that evolve on different timescale, which enables the separation of their dynamics into "slow" and "fast." This property is exploited by methods based on timescale separation to reduce the complexity of a system by approximating the fast variables with their steady-state values. This approach preserves the biological meaning of the state space variables and favors modularity between processes with widely separated timescales (Grunberg and Del Vecchio 2019). However, timescale separation methods are only applicable when the underlying biochemical reaction network behaves like a closed system and may otherwise yield reduced order models whose behavior is qualitatively different from that of the original system (Stoleriu et al. 2004(Stoleriu et al. , 2005Flach and Schnell 2006;Pedersen et al. 2008). A popular alternative to methods based on timescale separation is lumping (Wei et al. 1969;Tomlin et al. 1994;Koschorreck et al. 2007;Sunnåker et al. 2011;Rao et al. 2014), which aggregates state space variables into "macroscopic" variables to reduce the dimensionality of a system. While in some cases this maintains biological interpretability and compatibility with open systems modeling (Rao et al. 2014), lumping often requires expert knowledge, which makes it unappealing for constructive design and quantitative verification. Another common approach to model reduction of biochemical systems is based on sensitivity analysis and optimization (Danø et al. 2006;Zhang and Goutsias 2010;Prescott and Papachristodoulou 2012;Hangos 2013), which build reduced order models by minimizing an error function (such as the sensitivity of a variable) within a given range of candidate models. Despite its wide applicability, this approach can be highly demanding from a computational viewpoint for large-scale models and, in general, offers no a priori guarantees on the behavior of the reduced order model. A notable exception is given by methods based on balanced truncation (Liebermeister et al. 2005;Hardin and van Schuppen 2006;MeyerBase and Theis 2008;Sootla and Anderson 2014) and on the use of related Linear Matrix Inequalities (LMIs) (Anderson et al. 2011). In broad terms, balanced truncation first computes a change of coordinates in which the degrees of reachability and observability of each state are the same; the states which are least controllable and observable are then truncated to obtain a reduced order model. A nice feature of balanced truncation is that it preserves stability and that it offers a priori error bounds in the L 2 norm. However, balanced truncation is primarily a linear method and, hence, is not directly applicable to capturing important nonlinear behaviors observed in biology, such as multistability and limit cycle oscillations (Fall et al. 2005;Murray 2007;Alon 2007;Keener and Sneyd 2008;Del Vecchio and Murray 2015). Motivated and inspired by this line of research, the present paper extends the applicability of balanced truncation to behaviors that are not restricted to the stability of a single equilibrium.
The model reduction problem has been extensively studied also in the systems and control literature (Antoulas 2005). The problem is well understood for finite-dimensional, linear, time-invariant systems, for which standard methods are based on balanced truncation (Moore 1981;Glover 1984), on moment matching (Georgiou 1983;Kimura 1986;Antoulas et al. 1990;Georgiou 1999), and on a combination of both approaches (Gugercin and Antoulas 2006;Gugercin 2008;. Oscillatory behaviors can be approximated, in principle, using balanced truncation for linear periodic systems (Varga 2000) and linear time-varying systems (Sandberg and Rantzer 2004). However, both approaches lead to reduced models whose state space dimension may vary over time and which are thus not amenable to analysis and design. Several methods exist to approximate the local behavior of nonlinear, time-invariant systems around equilibrium (Antoulas 2005), including proper orthogonal decomposition (Berkooz et al. 1993), balanced truncation (for stable input-affine systems) (Scherpen 1993), empirical balanced truncation (Hahn and Edgar 2002), moment matching (Astolfi 2010; Padoan and Astolfi 2019) discrete empirical interpolation (Chaturantabut and Sorensen 2010), high-order moment matching (Asif et al. 2020) and H 2 -optimal model reduction (Benner et al. 2018). However, these methods do not provide a priori guarantees on the global behavior of reduced order models away from equilibrium, which limits their applicability to multistable and oscillatory systems. The problem of approximating the global behavior of a nonlinear system is indeed largely open. Besselink and co-authors have recently proposed a model reduction framework which preserves (incremental) stability properties for systems that can be decomposed as the feedback interconnection of a largescale stable linear system and a contractive nonlinear system (Besselink et al. 2009(Besselink et al. , 2013(Besselink et al. , 2014. In the recent papers (Padoan et al. 2021(Padoan et al. , 2020, we have generalized these results to multistable and oscillatory systems that can be decomposed as the feedback interconnection of a large-scale stable linear system and a dominant nonlinear system. The present paper extends our preliminary results to general dominant systems motivated by the increasing need of approximation tools for biological systems away from equilibrium.
Contributions This paper proposes a model reduction framework for biological systems whose behavior is not restricted to the stability of a single equilibrium. Motivated and inspired by the series of works (Besselink et al. 2009(Besselink et al. , 2013(Besselink et al. , 2014 and our preliminary results (Padoan et al. 2020(Padoan et al. , 2021, the present paper revisits the model reduction problem for nonlinear systems in light of differential analysis and dominance theory Miranda-Villatoro et al. 2018;Padoan et al. 2019aPadoan et al. , b, 2020Padoan et al. , 2021. A nonlinear enhancement of classical balanced truncation for linear stable systems is developed for p-dominant systems. Following the paradigm of classical balanced truncation, the model reduction problem is reduced to the simultaneous diagonalization of a pair of gramian matrices, computed using the linearization of a system and LMIs subject to a fixed inertia constraint. The asymptotic behavior of a reduced order model is then characterized by ensuring that the property of p-dominance is preserved. Although the biological meaning of state space variables is not necessarily preserved, our approach offers a series of benefits. First, it favors tractability leveraging on standard convex optimization tools to build reduced order models. Second, the quality of a reduced order model can be interpreted in quantitative terms using classical control-theoretic notions and tools, such as eigenvalues, Nyquist diagrams and Bode diagrams. Third, it provides a compositional model reduction framework, where the global emergent behavior of a complex biological system can be approximated using reduced order models of its elementary components. As a motivating example, we consider the problem of approximating the Goldbeter model for the expression of the per gene in Drosophila (Goldbeter 1995)-a paradigmatic example of biological oscillations (see Goldbeter 1996;Fall et al. 2005;Murray 2007;Keener and Sneyd 2008). Numerical results suggest that the proposed framework may be relevant to the approximation of a wide range of biological systems.
Paper organization The remainder of this work is organized as follows. Section 2 first recalls preliminary results on balanced truncation for stable linear systems and on dominance theory. Classical balanced truncation is then revisited in light of dominance theory to develop a model reduction method for the analysis of multistable and oscillatory systems. Section 3 illustrates the applicability of the proposed model reduction framework by means of a worked-out numerical example. Section 4 discusses main benefits and possible improvements of the proposed method. Section 5 concludes the paper with a summary of our main results and an outlook to future research directions. The Appendix contains more detail on the algorithms behind parameter selection for dominance analysis and simultaneous diagonalization of matrices with a fixed inertia.
Notation R and C denote the set of real numbers and the set of complex numbers, respectively. R + denotes the set of nonnegative real numbers. σ (A) denotes the spectrum of the matrix A ∈ R n×n . The inertia of the matrix P ∈ R n×n is defined as In(P) = (n − , n 0 , n + ), where n − is the number of eigenvalues of P in the open left half plane, n 0 is the number of eigenvalues of P on the imaginary axis and n + is the number of eigenvalues of P in the open right half plane, respectively.

Balanced truncation for linear stable systems
Consider a linear, time-invariant, system described by the equationṡ in which x ∈ R n , u ∈ R m , y ∈ R l , A ∈ R n×n , B ∈ R n×m , C ∈ R l×n , respectively. Assume that system (1) is stable and minimal, i.e, reachable and observable.

Balancing
Balancing for system (1) consists in finding a coordinates transformationx = T −1 x, with T ∈ R n×n and det(T ) = 0, such that the reachability gramian P ∈ R n×n and the observability gramian Q ∈ R n×n , defined implicitly by the Lyapunov equations are both diagonal and, if possible, equal. Stability of system (1) ensures existence and positive definiteness of the gramians, while minimality ensures that these are full rank (Antoulas 2005). A change of coordinatesx = T −1 x for system (1) acts on the reachability gramian and the observability gramian as Balancing thus amounts to finding a transformation T which simultaneously diagonalizes the positive definite matrices P and Q. This is a classical problem in linear algebra with a well-known solution (Bernstein 2009, p. 422).
in which case the corresponding representation in coordinates of system (1) is said to be (principal-axis) balanced (Antoulas 2005). The diagonal elements of are referred to as the Hankel singular values of the system. The Hankel singular values do not depend on the particular coordinate system and their magnitude quantifies the influence of each state on the overall input-output behavior of the system (Antoulas 2005).

Model reduction by balanced truncation
Balanced truncation for system (1) consists in eliminating, by truncation, the state variables corresponding to the least n − r Hankel singular values, where r ≤ n is the order of the reduced order model. The resulting reduced order model with ξ ∈ R r ,û ∈ R m ,ŷ ∈ R l , is guaranteed to be stable and to satisfy the L 2 error bound 1 for all u ∈ L 2 . While classical balanced truncation only applies to linear, time-invariant, systems, several nonlinear extensions exist (see, e.g., Scherpen 1993;Hahn and Edgar 2002;Besselink et al. 2009Besselink et al. , 2013Besselink et al. , 2014Benner et al. 2018). All those references however consider stable nonlinear systems with a unique equilibrium attractor, except for (Benner et al. 2018), which, however, does not provide a priori guarantees on the global behavior of reduced order models away from equilibrium. The goal of the present paper is to propose model reduction methods for nonlinear systems with multiple stable equilibria or stable limit cycle attractors. The proposed approach is based on recent results from dominance theory, which we discuss in the next section.

Dominance theory
Consider a nonlinear, time-invariant system and its variational dynamics described by the equationṡ in which x ∈ R n , u ∈ R m , y ∈ R l , f : R n → R n is a continuously differentiable vector field, B ∈ R n×m and C ∈ R l×n are constant matrices, δx ∈ R n , δu ∈ R m , δ y ∈ R l (identified with the respective tangent spaces), and ∂ f is the Jacobian of the vector field f .

Definition 1 [ p-dominance] (Forni and Sepulchre 2019)
The system (8) is said to be p-dominant with rate λ : R n → R + if there exist ε ∈ R + and a symmetric matrix P ∈ R n×n , with inertia In(P) = ( p, 0, n − p), such that the prolonged system (8)-(9) satisfies for all (x, δx) ∈ R n × R n and for δu = 0. The property is strict 2 if ε > 0.
The property of p-dominance yields the quadratic form which allows one to rewrite the inequality (10) aṡ The quadratic form V (δx) thus characterizes the contraction of the cones For further detail on dominance theory, the reader is referred to the paper (Forni and Sepulchre 2019).

Differential balancing
Differential balancing for system (8) consists in finding a change of coordinatesx = T −1 x, with T ∈ R n×n such that det T = 0, such that the reachability gramian P ∈ R n×n and the observability gramian Q ∈ R n×n , defined implicitly by the Lyapunov inequalities for x ∈ S ⊂ R n , are both diagonal and, if possible, equal.
Differential balancing is conceptually similar to classical balancing for linear, time-invariant, stable systems, but a few important points need to be highlighted. The Lyapunov inequalities (13) and (14) generalize the classical Lyapunov equations (2) and (3) and need to be solved uniformly in x in a given subset S of the state space. There are many ways to solve families of LMIs of the form (13) and (14). It is common practice to reduce families of LMIs like (13) and (14) to a finite family of LMIs through convex relaxation, see, e.g., (Boyd et al. 1994) and references therein. Another approach is to solve the inequalities (13) and (14) on a suitably dense grid covering the set S, which then allows one to exploit the continuity properties of the Jacobian matrix ∂ f (x) to ensure that the inequalities are satisfied everywhere. This approach is analogous to a classical computational procedure for finding a controllability Gramian and an observability Gramian for a parameter-dependent system (see, e.g., Wood et al. 1996 and Son and Stykel 2017 for recent developments). An indepth discussion on the computational aspects behind the solution families of LMIs of the form (13) and (14) is given in , Section VI). However, for the convenience of the reader, a heuristic procedure to select the parameters p and λ(x) is given in Appendix A.
In contrast to classical balancing, (13) and (14) do not necessarily admit a solution. Moreover, if solutions do exist, the reachability and observability gramians are not necessarily positive definite, but have a fixed inertia. The existence of solutions for (13) and (14) ensures p-dominance of the system (8), since both (13) and (14), together with B B T ≥ 0 and C T C ≥ 0, directly imply (10).
In analogy with classical balancing, a change of coordinatesx = T −1 x for system (8) naturally induces a change of coordinates δx = T −1 δx on its linearization (9). This, in turn, acts on the reachability gramian and the observability gramian as in (4). Differential balancing thus amounts to finding a transformation T which simultaneously diagonalizes the matrices P and Q. This problem admits a solution whenever the spectrum of P and Q satisfies certain assumptions (see, e.g., Uhlig 1973;Kenney and Hewer 1987;Therapos 1989). The solution of this problem and the corresponding algorithms are discussed in Appendix A.
Similar to classical balancing, the matrix T ∈ R n×n is said to be a (principal-axis) balancing transformation if (5) holds, in which case the corresponding representation in coordinates of system (8) is said to be (principal-axis) balanced. The diagonal elements of are said to be the characteristic values of the system.

Model reduction by differential balanced truncation
Differential balanced truncation for a p-dominant system (8) consists in eliminating, by truncation, the state variables corresponding to the least n −r characteristic values (in absolute value), where r ≤ p is the order of the reduced order model. To illustrate this method, assume that the Lyapunov inequalities (13) and (14) are solved by P ∈ R n×n and Q ∈ R n×n , with In(P) = In(Q) = ( p, 0, n − p). Further, assume that system (8) is (principal-axis) balanced, so that the gramians can be partitioned as with 1 = diag(σ 1 , . . . , σ r ) and 2 = diag(σ r +1 , . . . , σ n ) such that In( 1 ) = ( p, 0, r − p) and In( 2 ) = (0, 0, n − r ), respectively. Then (15) directly induces the partitions x = [ x T 1 x T 2 ] T , with x 1 ∈ R r and x 2 ∈ R n−r , and The reduced order model obtained by differential balanced truncation is defined by setting x 2 = 0 and discarding the dynamics of x 2 , yieldinġ with ξ ∈ R r ,û ∈ R m ,ŷ ∈ R l . The reduced order model (17) is (principal-axis) balanced and p-dominant with rateλ(ξ ) = λ(ξ, 0). To see this, consider the partition (16) and the Lyapunov inequalities (13) and (14) and shows that (17) is (principal-axis) balanced. Moreover, both (18) and (19) imply By Definition 1, the reduced order model (17) is thus

Results
We illustrate the proposed model reduction framework by approximating the behavior of a classical biological model of circadian oscillations: the Goldbeter model (Goldbeter 1995). Experiments have been conducted using standard routines of MATLAB (R2019b) on a 3.5 GHz Intel Core i7 processor.

The Goldbeter model
The Goldbeter model is a classical model in cellular physiology which describes circadian oscillations in the expression of the gene per in Drosophila (Goldbeter 1995) (see also Keener and Sneyd 2008, p. 440). The model is governed by the equationṡ in which M ∈ R + is the concentration of per mRNA, P 1 ∈ R + , P 2 ∈ R + and P 3 ∈ R + are the concentrations of unphosphorylated, monophosphorylated and biphosphorylated PER protein, and P N ∈ R + is the concentration of PER protein in the nucleus, respectively. With the original parameters-reported in Table 1 for the reader's convenience -the Goldbeter model has bounded solutions, a unique stable limit cycle and a unique unstable equilibrium point (Murray 2007).

Dominance analysis
To illustrate our model reduction framework, we consider system (21) with an added exogenous input u to (21a) and with P N as the output variable y, namely which can be described as a system of the form (8) by defining x = [ M P 0 P 1 P 2 P N ] T ∈ R 5 and System (22)  The dominance properties of system (22) have been analyzed selecting the parameters p and λ(x) by means of the heuristic procedure given in Appendix A, which, taking the limit cycle S as the region of interest, yields p = 2 and λ(x) = k 1 x 4 − k 2 x 5 − 0.4. Figure 1 shows the eigenvalues of the Jacobian of system (22)-given in (24)-at distinct samples {x(t k )} N k=1 of the solution at steady state. Note that ∂ f (x) has 2 eigenvalues with real part larger than −λ(x) (red) and 3 eigenvalues with real part less than −λ(x) (black) for λ(x) = k 1 x 4 − k 2 x 5 − 0.4, which suggests the existence of a dominant splitting. By evaluating the Jacobian of the system ∂ f (x) at distinct samples {x(t k )} N k=1 of the solution at steady state, the LMI (10) has been solved for ε = 0.001 and λ(x) = Re σ(∂f (x)) Im σ(∂f (x)) Fig. 1 Eigenvalues of the Jacobian of the system (22) at distinct samples {x(t k )} N k=1 of the solution at steady state. Note that ∂ f (x) has 2 eigenvalues with real part larger than −λ(x) (red) and 3 eigenvalues with real part less than −λ(x) (black) for λ( whose the inertia is In(P) = (2, 0, 3). The same analysis has been repeated for different sampling times obtaining similar results. In agreement with Theorem 1, this suggests that the Goldbeter model (22) is 2-dominant with rate λ(x) = k 1 x 4 − k 2 x 5 − 0.4 in a neighborhood of the limit cycle S.

Differential balanced truncation
We now turn to the question of building a reduced order model of system (22). First, the reachability and observability gramians have been computed by sampling the Jacobian of the system as above and solving for ε = 0.001 and λ(x) = k 1 x 4 − k 2 x 5 − 0.4 the LMIs (13) and (14) whose inertia is In(P) = In(Q) = (2, 0, 3).
The spectrum of the matrix P Q is real and positive: σ (P Q) = { 2847.3, 1590.1, 984.4, 161.3, 21.6 }. The system (22) thus admits a (principal-axis) balancing transformation by virtue of Theorem 2 (given in Appendix A). A (principal-axis) balancing transformation T ∈ R 5×5 has been computed using Algorithm 1 (also given in Appendix A), which yields The change of coordinatesx = T −1 x simultaneously diagonalizes the gramians according to (4), yieldinḡ For ε = 0.001 and λ(x) = k 1 x 4 − k 2 x 5 − 0.4, the characteristic values of the system are thus σ 1 = −39.877, σ 2 = −31.376, σ 3 = 4.651, σ 4 = 12.699, σ 5 = 53.360. Note that σ 3 and σ 4 are much smaller in absolute value relative to σ 1 , σ 2 and σ 5 , which suggests that the corresponding state variablesx 3 andx 4 have a negligible effect on the overall (differential) input-output behavior and, hence, may be eliminated. This intuition is confirmed by Fig. 2, which shows the time history of the solutions of system (22), with u = 0, in the new coordinatesx (top) as well as those of the reduced order models of order r = 4 (middle) and r = 3 (bottom) obtained using differential balanced truncation, respectively. Note that differential balanced truncation indeed eliminates the state variables which, in the new coordinates, have the least variation relative to the other state variables. Figure 3 (top) shows the time history of the output of system (22), with u = 0, in the new coordinates (solid) and of the reduced order models of order r = 3 (dotted) and r = 4 (dashed) obtained using differential balanced truncation, as well as those of the corresponding errors in absolute value (bottom), respectively. Note that while the output of the reduced order model of order r = 3 is out of phase with the oscillation of the original system, the output of the reduced order model of order r = 4 tracks well the output of the original system.  (22), with u = 0, (solid) and those of the reduced order models of order r = 4 (dotted) and r = 3 (dashed) obtained using differential balanced truncation, respectively. Bottom: Time history of the corresponding output errors in absolute value (logarithmic scale) Re σ(∂f 1 (ξ,0)) Im σ(∂f 1 (ξ,0)) Fig. 4 shows the eigenvalues of the Jacobian of the system (22), with u = 0, at distinct samples x(t k ) of the solution at steady state (top) and those of the corresponding reduced order models of order r = 4 (middle) and r = 3 (bottom) obtained by differential balanced truncation. The figure suggests that each reduced order model preserves 2-dominance with rate λ(x). It is interesting to note that differential balanced truncation eliminates the modes associated with the most negative eigenvalues of the linearized dynamics. This is consistent with the intuition that a "good" reduced order model should capture the dominant dynamics and, hence, eliminate the fastest transient modes of the (differential) input-output behavior of the system. Eigenvalues of the Jacobian of the system (22), with u = 0, at distinct samples x(t k ) of the solution at steady state (top) and those of the corresponding reduced order models of order r = 4 (middle) and r = 3 (bottom) obtained by differential balanced truncation, respectively

Discussion
The proposed model reduction framework offers a series of benefits. First, the construction of reduced order models relies on standard convex optimization and linear algebra tools, namely the solution of LMIs and the simultaneous diagonalization of a pair of matrices. This renders model reduction easy to implement and favors tractability. Second, the quality of a reduced order model can be interpreted in quantitative terms using classical control-theoretic notions and tools, such as eigenvalues, Nyquist diagrams and Bode diagrams. Third, our approach is compositional and compatible with open systems modeling. Combining our model reduction framework with the small-gain theorem for p-dominance , the global emergent behavior of a complex biological system can be approximated using reduced order models of its elementary components.
On the other hand, the proposed model reduction framework can be improved in several ways. A first limitation of our approach is that it relies on a change of coordinates which, in general, does not preserve the biological meaning of state space variables. This is a well-known drawback of model reduction methods based on balanced truncation (Snowden et al. 2017), which can be circumvented by imposing a given sparsity pattern while computing the reachability and observability gramians as discussed, e.g., in Sootla and Anderson (2017). A similar approach can be taken within our framework to maintain the biological interpretation of state space variables by requiring that the Lyapunov inequalities (13) and (14) have a given sparsity pattern.
A second limitation is that the LMIs (13) and (14) are only verified in a neighborhood of the limit cycle, resulting in a local result. While our analysis is not infinitesimal and can be in principle adapted to any desired region of the state space, it would be of interest to compare the proposed results to the model reduction of a periodic linear system obtained by linearizing the nonlinear dynamics along the limit cycle. Model reduction for periodic linear systems has been addressed, e.g., in Varga (2000); Sandberg and Rantzer (2004).
A third limitation of our approach is that it only applies to systems of moderate size. Indeed, the dimension of a system directly influences the number of variables and LMIs required to solve (13) and (14), which is limited to a few thousands in currently available solvers. Nonetheless, the compositional nature of our framework enables one to partition a large-scale system into smaller subsystems, compute reduced order models for each individual subsystem and finally obtain by interconnection an overall reduced order model, the behavior of which can be characterized a posteriori using the small gain for p-dominance .
A further limitation concerns the computational cost of the simulation of the reduced order model (17). Given a balancing transformation T ∈ R n×n , the simulation of the reduced order model (17) requires the evaluation of f 1 (ξ, 0) = S 1 f (S 1 ξ), where S 1 are the first r columns of T andS 1 are the first r rows of T −1 , respectively. This means that, in general, the computational cost of evaluating f 1 (ξ, 0) depends on the order n of the original system. As a result, the reduced order model may offer a limited performance speed-up compared to the original model. This computational issue has been previously addressed in the literature (see, e.g., Chaturantabut and Sorensen 2010). Nonlinear functions can be projected onto a subspace that approximates the space generated by the nonlinear terms and that is spanned by a basis of lower dimension (see Chaturantabut and Sorensen 2010, Section 2.2 for further detail). The same strategy can be used to reduce the computational cost of the simulation of the reduced order model (17), since our framework also relies on constant projections.
Finally, our framework does not provide an a priori error bound. Similar to the case of linear time-invariant stable systems, it is reasonable to expect that the approximation error is bounded by a function of the neglected characteristic values. A thorough analysis of this issue is beyond the scope of this preliminary study and is the subject of ongoing research. Nonetheless, we observe that in the context of our example the differential input-output behavior of the original system is well approximated. This can be appreciated by considering the Nyquist diagrams of the family of transfer functions defined by the linearization of the original system as well as those of the family of transfer functions defined by the linearization of each reduced order model which are formed by tracing s ∈ C around the Nyquist "D contour" consisting of the imaginary axis combined with an arc at infinity connecting the endpoints of the imaginary axis (see Åström and Murray 2020 for further detail). Nyquist diagrams are widely employed in control theory (Åström and Murray 2020) and provide considerable insight into the input-output behavior of a dominant system (Miranda-Villatoro et al. 2018;Padoan et al. 2019a, b). Figure 5 shows the Nyquist diagrams of the family of transfer functions (30) and (31) defined by the original system (top) and by the reduced order models of order r = 4 (middle) and r = 3 (bottom), respectively. Not only the Nyquist diagrams of each member of the family of transfer functions (30) and (31) have a similar shape, but the distance between one another is negligible. The quality of the approximation is reflected in Fig. 6, which shows the Bode diagrams the magnitude of the family of error transfer functions defined as Note that the approximation the error is bounded as for the reduced order model of order r = 4 and as for the reduced order model of order r = 3, which indicates that the differential input-output behavior of the original  (30) and (31) associated with the linearizations of the original system (top) and the reduced order models of order r = 4 (middle) and r = 3 (bottom) system is well approximated by that of both reduced order models.

Conclusion
The paper has outlined a model reduction framework for biological systems whose behavior is not restricted to the stability of a single equilibrium. Classical balanced truncation for linear stable systems has been extended to dominant nonlinear systems using differential analysis. The asymptotic behavior of reduced order models has been characterized by ensuring that the property of p-dominance is preserved. This approach is tractable and offers quantitative tools to predict the behavior of a reduced order model a priori. Preliminary numerical results suggest that the proposed model reduction framework may be relevant to the approximation of a wide range of biological systems.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Algorithms
This section provides a heuristic procedure for the selection of the parameters p and λ(x) in dominance analysis as well as necessary and sufficient conditions for the existence of a (principal-axis) balancing transformation and an algorithm to find such a transformation.

A.1 Selecting p and (x) for dominance analysis
A necessary condition for p-dominance of the system (8) is that the spectrum of the family of matrices ∂ f (x) + λ(x)I admits a uniform splitting into p eigenvalues with positive real part and n − p eigenvalues with negative real part (Forni and Sepulchre 2019, Section VI). The analysis of the spectrum of the Jacobian matrix ∂ f (x) can therefore inform the selection of the parameters p and λ(x). In particular, the following heuristic procedure can been used to estimate p and λ(x) based on the evaluation of the Jacobian matrix ∂ f (x) over a sufficiently fine grid of points {x k } N k=1 , with x k ∈ S and a subset S ⊂ R n of interest.

For each
The heuristic procedure outlined above comprises four basic steps. The first step consists in computing the spectrum of the Jacobian matrix ∂ f (x k ) at each point x k . The second step consists in estimating the parameter p based on the existence of a splitting into p "dominant" eigenvalues λ k,1 , . . . , λ k, p and n − p "transient" eigenvalues λ k, p+1 , . . . , λ k,n . Note that if such splitting does not exist, the system is not p-dominant for any p by the necessary condition mentioned above. The third step definesλ k as the negative part of the midpoint between Rλ k, p (the real part of the smallest dominant eigenvalue) and R λ k, p+1 (the real part of the largest transient eigenvalue). Finally, an estimate of the function λ(x) is computed by fitting {λ k } N k=1 to the samples {x k } N k=1 .

A.2 Balancing gramians with inertia
Recall that the goal of differential balancing for system (8) boils down to finding a transformation T ∈ R n×n which simultaneously diagonalizes the reachability gramian P ∈ R n×n and the observability gramian Q ∈ R n×n , defined implicitly by the Lyapunov inequalities (13) and (14). The following result provides necessary and sufficient conditions for the existence of a (principal-axis) balancing transformation T . The proofs are omitted as the claim follows directly from (Uhlig 1973), Corollary 1.4 see also (Kenney and Hewer (1987), Theorems 1 and 2) and (Therapos 1989). (13) and (14) admit solutions P ∈ R n×n and Q ∈ R n×n , with In(P) = In(Q) = ( p, 0, n − p). Then there exists a (principal-axis) balancing transformation if and only if the product P Q is similar to a (positive) real diagonal matrix.
We now provide a procedure to compute a (principalaxis) balancing transformation T ∈ R n×n , described in pseudo-code in Algorithm 1. The procedure comprises four basic steps. The first step (line 1) computes the eigenvalue decomposition S −1 P QS = 2 of the product P Q. The second step (line 2) defines the matrices X = S T QS and Y = S −1 P S −T , which can be simultaneously diagonalized since XY = Y X. The third step (line 3) computes a (unitary) coordinate transformation U which simultaneously diagonalizes X and Y , yielding X = U o U T and Y = U r U T . Finally, a (principal-axis) balancing transformation is computed as T = SU ( r −1 o ) 1/4 (line 4). The following result can be proved by direct computation.
Theorem 3 Suppose (13) and (14) admit solutions P ∈ R n×n and Q ∈ R n×n , with In(P) = In(Q) = ( p, 0, n − p), such that the product P Q is similar to a (positive) real diagonal matrix. Then the output of Algorithm 1 is a (principal-axis) balancing transformation.

Algorithm 1
Input: The gramians P ∈ R n×n and Q ∈ R n×n . Output: A (principal-axis) balancing transformation T ∈ R n×n . Assume: The solutions P ∈ R n×n and Q ∈ R n×n of (13) and (14) are full rank and such that the product P Q is similar to a (positive) real diagonal matrix. 1: Compute a coordinates transformation {S ∈ R n×n } such that S −1 P QS = 2 , with ∈ R n×n a (positive) real diagonal matrix. 2: Define X = S T QS and Y = S −1 P S −T . 3: Compute a (unitary) coordinates transformation U ∈ R n×n such that X = U o U T and Y = U r U T . 4: return T = SU ( r −1 o ) 1/4 .