A combined model reduction algorithm for controlled biochemical systems

Background Systems Biology continues to produce increasingly large models of complex biochemical reaction networks. In applications requiring, for example, parameter estimation, the use of agent-based modelling approaches, or real-time simulation, this growing model complexity can present a significant hurdle. Often, however, not all portions of a model are of equal interest in a given setting. In such situations methods of model reduction offer one possible approach for addressing the issue of complexity by seeking to eliminate those portions of a pathway that can be shown to have the least effect upon the properties of interest. Methods In this paper a model reduction algorithm bringing together the complementary aspects of proper lumping and empirical balanced truncation is presented. Additional contributions include the development of a criterion for the selection of state-variable elimination via conservation analysis and use of an ‘averaged’ lumping inverse. This combined algorithm is highly automatable and of particular applicability in the context of ‘controlled’ biochemical networks. Results The algorithm is demonstrated here via application to two examples; an 11 dimensional model of bacterial chemotaxis in Escherichia coli and a 99 dimensional model of extracellular regulatory kinase activation (ERK) mediated via the epidermal growth factor (EGF) and nerve growth factor (NGF) receptor pathways. In the case of the chemotaxis model the algorithm was able to reduce the model to 2 state-variables producing a maximal relative error between the dynamics of the original and reduced models of only 2.8% whilst yielding a 26 fold speed up in simulation time. For the ERK activation model the algorithm was able to reduce the system to 7 state-variables, incurring a maximal relative error of 4.8%, and producing an approximately 10 fold speed up in the rate of simulation. Indices of controllability and observability are additionally developed and demonstrated throughout the paper. These provide insight into the relative importance of individual reactants in mediating a biochemical system’s input-output response even for highly complex networks. Conclusions Through application, this paper demonstrates that combined model reduction methods can produce a significant simplification of complex Systems Biology models whilst retaining a high degree of predictive accuracy. In particular, it is shown that by combining the methods of proper lumping and empirical balanced truncation it is often possible to produce more accurate reductions than can be obtained by the use of either method in isolation. Electronic supplementary material The online version of this article (doi:10.1186/s12918-017-0397-1) contains supplementary material, which is available to authorized users.

1 Overview of model reduction methods 1

.1 The Galerkin projection
Methods of model reduction can be considered as a projection of the state-variables to a lower dimensional subspace V : dim (V) =ñ of the original phase-space, within which some relevant set of the system's trajectories can be adequately approximated. Mathematically, applying such a projection to obtain a reduced dynamical system is underpinned by the Galerkin projection [1] which will be introduced here.
Consider a basis B of the subspace V such that B = [b 1 , . . . , bñ] ∈ R n×ñ . Assuming B has been selected such that it provides an adequately accurate approximation of the original states x(t) within the subspace V, then withx(t) ∈ Rñ representing the reduced set of state-variables. Substituting this approximation into the stoichiometric model form yields where it is assumed that B is time-invariant. Additionally, ρ(t) ∈ R n is termed the residual and addresses the discrepancy emerging from the fact that Bx is typically not an exact solution of the system for all times. Now let W represent a subspace that is orthogonal to the residual ρ(t) with a basis C ∈ R n×ñ such that C ρ(t) = 0. Hence, left multiplying equation (2) by C produces Assuming C B is non-singular, this finally leads to a reduced dynamical system of the form This simplification of a dynamical system to a lower dimensional subspace is known as the Petrov-Galerkin projection. In the special case where B = C it is known simply as the Galerkin projection. In that case Such thatB is a generalised left inverse of B andBB = Iñ (theñ dimensional identity matrix). Whilst the explanation given above provides an account of how to apply a Petrov-Galerkin projection, it does not provide a methodology for finding suitable bases B and C for a given model.

Conservation analysis
To understand the nature of conservation relations and how they might be found computationally, first note that a common means for representing the network structure underlying a system of chemical equations is that of the stoichiometry matrix. This is an n × m matrix S, with each of the rows corresponding to a single species and each of the columns to a reaction. The matrix is populated such that its entries s ij give the net value of the stoichiometric coefficients (product minus reactant) of the i-th species in the j-th reaction. If the concentration of a particular species is not affected by a reaction the corresponding entry is populated with a 0. Hence, the sign of the entry indicates whether the species is a net reactant or a net product in the relevant reaction. A positive sign implies that the species is a product (i.e. the number of molecules is increased by the reaction), whilst a negative sign indicates that the species is a reactant (i.e. the number of molecules is decreased).
This matrix can be considered as mapping the vector of reaction rates v(x(t)) to the change in species concentration. Hence it is possible to represent the system of ODEs in the forṁ Now turning to conservation analysis and following the ideas outlined by Reder [2], the existence of conservation relations in a model implies that where Γ is an h × n matrix that will be referred to here as the conservation matrix, the rows of which represent the linear combinations of species that are constant in time. Alternatively, by integration, the h individual elements of which are known as conservation relations, with c ∈ R h representing a set of constants known as conserved values.
It is hence possible to express a model containing conservation relations in the form of a system of differential algebraic equations (DAEs). To see this, first partition x into two subsets: x d an h dimensional subset of the species with each element corresponding to a single species involved in a given conservation relation termed the dependent species. And x i an n − h dimensional subset accommodating all remaining state-variables, termed the independent species, such that Then from equation (8) Γ This is a system of linear equations and hence if Γ is expressed in reduced row echelon form, such that with I h representing the h dimensional identity matrix and N 0 a h × (n − h) matrix, it becomes apparent that This implies that the subset of dependent species x d can be eliminated from the governing system of ODEs by substituting in the appropriate element of equation (12). Hence, given the stoichiometric form given in equation (6), a system exhibiting conservation relations can be expressed in the form of a semi-explicit system of DAEs, such thatẋ where equation (13b) has been exploited in equation (13a) to obtain a system of ODEs such that statevariables x d are no longer explicitly given. Additionally, S i here represents the rows of the stoichiometric matrix corresponding to the independent state-variables x i . Obtaining the conservation matrix Γ, particularly for large systems, is often not feasible from simple inspection. To understand a more algorithmic approach for obtaining this matrix, begin by recalling the stoichiometric form of a model. Decomposing the stoichiometric matrix via the same partition as the set of species leads to the system ẋ d (t) However, via differentiation of equation (12) Hence, S d = −N 0 S i and therefore each conservation relation can be seen as corresponding to a linear dependency in the stoichiometry matrix. As such, conservation relations can be found by seeking the left null space Z n of S (i.e. via finding the null space of S T ) such that and hence Z n S = 0. This implies that and therefore via comparison to equation (7) it is clear that such that the conservation matrix is equal to the transpose of the left null space of the stoichiometry matrix. Hence, a more mathematically rigorous approach for finding conservation relations has been provided. For very large systems (i.e. n > 100) S will be a large and, typically, sparse matrix. As a result, solving the system of linear equations S z = 0 for each conservation relation may not always be numerically stable or efficient under traditional approaches such as Gaussian elimination. Therefore, the combined algorithm employs QR factorisation via Householder reflections, as is discussed in the main text.

Empirical balanced truncation
The description of empirical balanced truncation provided here is based upon the procedure given by Hahn and Edgar [3]. The aim is to construct two covariance matrices, known as the empirical controllability and observability Gramians, via repeated simulations of the system under perturbations of the input and the initial conditions respectively. As in the linear case, the empirical controllability Gramian provides information on how changes in the input will alter the state of the system. The empirical observability Gramian provides information on the total magnitude of the output that any given initial condition of the system can produce.
To construct the empirical Gramians, first recall that n refers to the original number of state-variables and l to the number of inputs. Additionally, we here define s to refer to the number of perturbation magnitudes and r to refer to the number of perturbation directions. Given this, we now define the following terms: • E n C = {e C,1 , . . . , e C,l } representing the set of standard unit vectors in the input space R l ; . , e O,n } representing the set of standard unit vectors in the state space R n ; • P r C = {P C,1 , . . . , P C,r }, a set of r orthogonal (P C,i P C,i = I, ∀i) l × l matrices representing the perturbation directions explored for the set of inputs; The Gramians are then computed using finite trajectory samples measured at regularly placed, discrete time points separated by a sampling interval ∆t. Additionally q is defined to be the number of points in the sample. To ensure that accurate empirical Gramians are obtained it is required that the system is approximately at equilibrium at some time t * ≤ q∆t. Using these discrete time point samples the controllability Gramian is defined as where ss and x ihm (t) is the state or condition of system at time t corresponding to the perturbed input u(t) = c m P C,h e C,i δ(t) + u 0 . Additionally, u 0 is the steady-state or unperturbed level of the input and x ihm ss represents the steady-state of the system produced under the perturbed input.
Meanwhile, the observability Gramian is similarly defined as where Φ hm ij = y ihm (t) − y ihm ss y jhm (t) − y jhm ss and y ihm (t) is the output of the system corresponding to the perturbed initial condition x ihm (0) = c m P O,h e O,i + x 0 , with x 0 representing the original, unperturbed initial condition of the system. Additionally, y ihm ss represents the steady-state of the output under this perturbed initial condition.
Once the Gramians have been computed the aim is to construct a balancing transformation of the statevariables which also acts to equalise and diagonalise the Gramians. This can be achieved by following the approach outlined in Methods Section of the main text.

Chemotaxis model summary
The specific model of chemotactic signalling analysed in the main text was originally described in Tindall et al. [4]. It is an 11 dimensional, nonlinear model describing a system of 12 biochemical reactions given by CheA + CheY k9 k−9 CheA · CheY, (22d) Whilst a more detailed account of this network can be found in the original paper, the model can be broadly understood as follows: Here, CheA represents a histidine kinase whose rate of autophosphorylation is modulated via extracellular attractant-receptor binding (in this case represented in equations (22a) and (22c) by the reaction rates k 1 (L) and k 10 (L), where L represents the concentration of extracellular ligand). An increase in attractant binding results in a decreased rate of phosphorylation of CheA. Two response regulator proteins (here represented by CheB and CheY and their interactions in equations (22b), (22d), and (22h)) compete to bind with protein CheA. Once bound CheA P will transfer its phosphoryl group, hence providing the only means of phosphorylation for proteins CheB and CheY. In the case of protein CheY this phosphotransfer is also reversible via the reactions given in equation (22c). After they have been phosphorylated both response regulator proteins steadily auto-dephosphorylate as given by equations (22f) and (22g). In the case of protein CheY, however, this process can be greatly accelerated by forming a complex with the phosphatase CheZ as described by equation (22e). Note that it is the concentration of phosphorylated protein CheY that regulates the process of chemotaxis by binding with parts of the flagellar motor complex (the process is not explicitly described in this model). The model differentiates itself somewhat from previously published approaches in that explicitly describes the intermediary complexes involved in the phosphorylation of proteins CheB and CheY.
Employing the Law of Mass Action, chemical equations (22) can be modelled by the following system of ODEs , and x 11 (t) = [CheB P ]. The initial condition employed in analysing this system was The specific parameterisation employed for this model can be found in Table 1.
To understand the input u, note that receptor binding of the extracellular chemotactic ligand (denoted by L) is described by Michealis-Menten kinetics, such that where K is a Michaelis-Menten constant and h is a Hill-coefficient. For the sake of simplicity the input u(t) is thus defined here to be such that k 1 (L) = k 10 (L) = k 1 u(t).  [4].

Reduction
In the following sections we will describe how the various steps of the combined model reduction algorithm were applied to the model of bacterial chemotaxis and reproduce the results it obtained.

.1.3 Lumping
The algorithm then proceeds to reduce the system via proper lumping to a 5 dimensional form. In this section we will symbolically reproduce this reduced model. Through the forward selection procedure described in the main text, the algorithm determines that the optimal scheme is to lump together state-variables z 3 (τ ), z 4 (τ ), and z 8 (τ ) as one. These variables correspond to the original species CheA P · CheY, CheA · CheY P , and CheY P , respectively. If we assume a vector of state-variables, such that y(τ ) = [y 2 (τ ), y 3 (τ ), y 4 (τ ), y 6 (τ ), y 7 (τ ), y 8 (τ ), y 10 (τ )] , then this scheme corresponds to the lumping matrix Hence we have a set of reduced state-variablesỹ(τ ) = Ly(τ ) such that Following the general descriptions provided in the main text it is now necessary to construct a generalised right inverseL of the lumping matrix prior to application of the Galerkin projection. Such a matrix can be constructed by computingL where X represents a diagonal matrix, typically containing steady-state values of the system. In the case of the nondimensionalised model described above, this would yield a matrix of the form Given these matrices it is then possible to construct a reduced description of the system's dynamical where the new parameters (γ 1 , γ 2 , and γ 3 ) are defined in Table 4. The specific values reported were calculated via the novel, averaged inverse approach that is outlined in the main manuscript.

Empirical Balanced Truncation
Given the system reduced via the application of nondimensionalisation, conservation and lumping described by equations (33) and the parameter sets defined by Tables 2, 3 and 4, it is now possible to apply empirical balanced truncation in order to achieve further reduction of this system following the combined model reduction algorithm.
To construct the empirical Gramians we employed simulated data for 100 distinct perturbations to both the unpertubed model input u 0 = 1 and each of the unperturbed system's initial conditionsỹ 0 , in all instances perturbations were uniformly sampled from of 0.2 to 1.8 times the parameter in question's original or unperturbed values. Following the approach described in Section 1.3, this yielded the Gramians Calculating the balancing transformation following the procedure given in the main text then yields If we now seek to create the 2 dimensional reduced system described in the main text we first need to define reduced variablesx(τ ) = [x 1 (τ ),x 2 (τ )] . However, to increase the accuracy of this reduced model it should be residualised, such that for reconstructing the values of the set of nondimensionalised variables y(τ ) from this reduced model we have Given these state-variables, the reduced model is then with associated initial conditionsx 1 (0) = 1.6775 andx 2 (0) = 1.4560.

ERK activation model
Due to the size of the ERK activation model it is not possible to provide the same level of detail into its reduction as was provided for the model of bacterial chemotaxis signalling. We have, however, provided a number of Matlab files that describe the form of the model at various stages of reduction along with its initial conditions, associated lumping matrices and their generalised inverses, the empirical gramians and the balancing transformation. Use of these files does require Matlab's Symbolic Math toolbox. All of this is made available in Additional file 2 -ERK Activation Reduction Files.

Combined algorithm overview
In this section a more detailed account of the combined model reduction algorithm is provided. Figure 1 depicts the high level steps of the overall algorithm. The following sections correspond to each of these steps and describe the specific process and implementation in Matlab by which they are achieved. This section provides a more detailed unpacking of the algorithm steps given in Figure 2 of the main text. Figure 1. A shortened overall schematic of the combined model reduction algorithm. This is a submodule, schematic representation of the combined model reduction algorithm that will be detailed in sections 3.1 to 3.5.

Importing SBML
• Have the user define an SBML file representing the model, referred to here as 'SBML.xml'. Additionally, define the desired input values u and outputs y(t) in terms of the species given in the model. Finally, define a 'natural' or unperturbed input level u 0 which corresponds to the initial condition of Figure 2. Import and symbolic conversion process in the combined model reduction algorithm.
the pathway before any intervention occurs (for example the concentration of an extracellular ligand considered as an input might be negligible until the administration of the particular drug of interest, hence u 0 = 0).
• The algorithm uses the LibSBML Matlab library along with the open-source SBToolbox2 package to import SBML.xml as an SBmodel object named 'mod'.
• Commands from the SBToolbox2 package are then additionally used to extract specific model information from 'mod' as variables into the Matlab workspace. Specifically: the set of species named 'Species'; the set of parameters named 'Params'; the set of reactions named 'Reacts'; and the stoichiometry matrix named 'N '.
• The model is simulated with the initial conditions given in the SBML file and under the undisturbed input u 0 . Once steady-state has been reached this undisturbed equillibrium is defined to be the new initial condition of the system.
• A symbolic representation of the model named 'DiffFuncs' and its variables is then constructed via the Symbolic Matlab package. To achieve this the algorithm defines a set of symbolic variables x i (t) of equal length to the vector Species; creates a symbolic vector of reactions 'SReacts' by substituting Species(i) = x i (t) into Reacts; -Redefines the outputs y(t) in terms of these new state-variables x i (t) via substitution; and computes DiffFuncs as the product of the stoichiometry matrix and SReacts, such that DiffFuncs = N · SReacts.
• DiffFuncs can then be converted to a form for simulation via the functions 'odeToVectorField' and 'matlabFunction'.
• For future reference, now simulate the system under all desired input values u and calculate the outputs y(t) under each. Store these results as the original output of the system 'OutO'.

Nondimensionalisation
• Begin by constructing a complete list of the parameters and their associated units: construct a list of each monomial term contained within SReact named SR. Here the inbuilt 'strsplit' function is employed to create the list of monomials and x i (t) = 1 is used for all i so as to access only the coefficient of each monomial; Figure 3. Nondimesionalisation process in the combined model reduction algorithm.
additionally, compute the degree of each monomial within SR; use this list of degrees to inform the fundamental units associated with the kinetic rate coefficient of each monomial. Here use T to represent temporal units and M to represent mass based units; and if the conserved totals are not given in the predefined model parameters Params, then it is necessary to calculate them here via the steps outlined in Section 1.1.
• Substitute scaled variables, such that x i (t) = α i z i (t) and t = βτ , into the symbolic system of ODEs (DiffFuncs).
• Compute a list of the new coefficients of each monomial in this scaled system of ODEs (NDP).
• Use the list of units to create a list of possible nondimensionalisations for the system (the algorithm is restrained to [β] = T such that β is the reciprocal of a 1/T kinetic rate coefficient and [α i ] = M, ∀i such that each α i is chosen from the list of conserved total masses. Other nondimensionalisations for rate terms associated with nonlinear monomials will exist, but the algorithm is currently limited in this regard for the sake of computational efficiency).
• Even for modestly sized systems (∼10's of ODEs) this will produce a very large list of possible nondimensionalisations. Hence, instead of evaluating them all, randomly select a sample of 50 possible nondimensionalisations (or use the total number if this is less than 50).
• Loop through this possible list of rescalings, substituting each into NDP with individual rescaling referred to as NDP(i).
• Substitute this nondimensionalisation into the symbolic system of ODEs to be regarded as the new model.

Conservation analysis & application
• Use an empirical Jacobian estimation algorithm around the initial conditions of the model. Here this is achieved using the SBtoolbox2 command 'SBjacobian'.
• Use this Jacobian estimate to approximate the speed of each reactant by simply using the diagonal elements. Rank the species in terms of these elements with larger values corresponding to faster species.
• Apply Matlab's inbuilt QR factorisation and reduced row echelon form algorithms to obtain the conservation relations for the model. Essentially this follows the approaches outlined in Section 1.2; apply QR factorisation 'qr' with pivoting to the stoichiometry matrix N to yield an orthogonal matrix Q, an upper triangular matrix R, and a pivot matrix P . Use the 'rref' function to place R in reduced row echelon form, and finally re-transform this matrix under P to yield the previously defined conservation matrix Γ.
• Substitute in initial conditions to obtain values for conserved totals.
• Solve conservation relations for particular species, favouring those species approximated to be the slowest.
• Substitute solutions into DiffFuncs to give a reduced set of ODEs (dismissing those differential equations corresponding to state-variables that have been reduced to algebraic relations). This reduced set of ODEs is now treated as the new model. • Begin by setting the 'previous lumping matrix' (L Prev ) and its generalised inverse (L Prev ) equal to the n dimensional identity matrix I n .

Proper lumping routine
• Create table to store lumping results, named here as 'ResTab'.
• Enter the first loop which is designed at each step to reduce the model by one further dimension.
Here, iterate for i from 1 to (n − 1).
-Calculate the new reduced dimensionñ = n − i and set the 'previous dimension' variable n prev = n + 1.
-Create a list of reduced symbolic variables z i (t), for i = 1 . . .ñ.
-Load previous lumping matrix (L Prev ) and its generalised inverse (L Prev ).
-Create a list of all possible proper, linear lumpings to dimensionñ based upon a one dimensional lumping of L Prev (there will be 1 2 (n prev [n prev − 1]) possible lumpings to trial as it represents all possible pairings of the current n prev state-variables). This list is referred to here as 'TempTab'.
-Enter a parallelised second loop using Matlab's 'parfor' function. Here, each thread trials a possible lumping from TempTab and calculates the error and stiffness coefficient associated with it. Iterate for j from 1 to size(TempTab). * The one dimensional n Prev toñ lumping matrix is defined by the j-th entry of TempTab, and is constructed here as L Temp . * Then the new overall lumping matrix is computed as L = L Prev · L Temp . * To trial both the steady-state and average lumping inverses,L andL respectively, now enter a final loop. Here, iterate for k from 1 to 2. · If k = 1: Define X = diag L Prev x * u0 , such that X represents an n prev × n prev diagonal matrix of the steady-state values x * u0 of the system under the unperturbed input.
· CalculateL Temp = XL Temp L Temp XL Temp · If k = 2: Use the novel average inverse defined in main text, such thatL Temp =L. · Now compute the overall lumping inverse such thatL =L Temp ·L Prev . · Construct the reduced system asż(t) = Lf (Lz(t)) by substituting 'reconstructed' state-variablesLz into DiffFuncs and then lumping this vector of symbolic equations via L. · Simulate this reduced model and calculate error as compared with the output of the original system (OutO). Additionally, re-estimate Jacobian and hence the stiffness coefficient χñ. Record both the error and stiffness for the specific lumping in TempTab. * End the lumping inverses loop.
-Select the lumping regime within a small error bound of the minimum error with the lowest stiffness coefficient.
-Set L Prev = L andL Prev =L.
-Record these 'optimal' error and stiffness results for reductionñ in ResTab.
• End main loop.
• Select the lowest error lumping below the threshold stiffness coefficient (set currently to χñ < 250 based upon numerical experimentation) and save the corresponding reduced model as the new system.
• Let n L represent the new dimensionality of the model after lumping.

Empirical balanced truncation routine
Compute the empirical controllability Gramian P following the approach outlined in Section 1.3: • Define P as an empty n L × n L matrix.
• Define a list of s random input perturbation sizes (M ). Here, typically 100 random sizes from 0 to the magnitude of u are selected.
• Define a list of r perturbation directions in input space. Here the algorithm currently employs r = 2 corresponding to an increase in all inputs and a decrease in all inputs. Additionally, define P (h) to represent an l × l diagonal matrix of the h-th perturbation direction.
• Define a set of l unit vectors in input space (E). Recall that l simply represents the number of inputs in the model.
• Enter the first loop going through each of the perturbation sizes. This loop is parallelized using Matlab's 'parfor' function. Here, iterate through m from 1 to s.
-Enter the second loop going through each of the perturbation directions. Here, iterate through h from 1 to r. * Enter the third loop going through each of the inputs. Here, iterate through i from 1 to l. · Define a perturbed input u p = M (m)P (h)E(i) + u ss . · Simulate the system under the perturbed input u p at q time-points equally spaced by a difference of ∆t, to yield x ihm (t). · Additionally define the steady-state of this simulation as x ihm ss . · Calculate P Temp = ∆t x ihm (t) − x ihm ss x ihm − x ihm ss . · Finally, add the contribution to the overall Gramian as P = P + 1 rsM (m) 2 P Temp * End the third loop.
-End the second loop.
• End the first loop.
• Return the empirical controllability Gramian P. Now compute the empirical observability Gramian Q following the approach outlined in Section 1.3: • Define Q as an empty n L × n L matrix.
• Define a list of s random input perturbation sizes (M ). Here, typically 100 random sizes from 0 to the magnitude of x 0 are selected.
• Define a list of r perturbation directions in state-space. Here the algorithm currently employs r = 2 corresponding to an increase in all initial conditions and a decrease in all initial conditions. Additionally, define P (h) to represent an n L × n L diagonal matrix of the h-th perturbation direction.
• Define a set of n L unit vectors (E) corresponding to each of the state-variables.
• Enter the first loop going through each of the perturbation sizes. This loop is parallelized using Matlab's 'parfor' function. Here, iterate through m from 1 to s.
-Enter the second loop going through each of the perturbation directions. Here, iterate through h from 1 to r. * Define a temporary matrix W for storing results from the following loop. * Enter the third loop going through each of the inputs. Here, iterate through i from 1 toñ. · Define a perturbed initial condition x p = M (m)P (h)E(i) + x 0 . · Additionally, recalculate the conserved totals using the perturbed initial conditions (so as to avoid violating conservation). · Simulate the system under the undisturbed input u ss with the perturbed initial condition x p at q time-points equally spaced by a difference of ∆t, to yield an output y ihm (t). · Additionally define the steady-state of the output of this simulation as y ihm ss . · Calculate w = y ihm (t) − y ihm ss . · Concatenate each column of w into the i-th row of W . * End the third loop. * Calculate Q Temp = W W . * Finally, add the contribution to the overall Gramian as Q = Q + ∆t rsM (m) 2 P (h)Q Temp P (h) .
-End the second loop.
• End the first loop.
• Return the empirical observability Gramian Q.
Given P and Q, now compute the balancing transformation as follows: • Perform a Cholesky factorisation of both of the Gramians using Matlab's inbuilt 'chol' function. This yields P = L L and Q = R R.
• Finally this yields a balancing transformation T = Σ − 1 2 V R and its inverseT = L U Σ − 1 2 . Now seek to truncate the system one dimension at a time as follows: • Define new reduced dimensionñ < n L .
• Calculate the transformed initial conditionsx 0 = T x 0 , such thatx 0,i ∈x 0 represents the i-th transformed initial condition.
• Simulate this reduced system and compare to the original output (OutO) to obtain a measure of error.
Finally, return the lowest dimensional system that maintains an acceptable maximal relative error bound ε < ε c . In the paper this has been set such that ε c = 0.05.