## 2 Materials and Methods
### The New Approach to Individual Subsystems Models: General Aspects
The second step consists of finding a minimum feedback arc set for each group (sources and sinks not included) which determines a set of iteration variables for each group of equations. Finally, if no partitioning is possible (i.e., there is only one group coincident with the whole ISMs), the node ordering procedure outlined below is applied directly to the RISM. In regards to nodes ordering, [PERSON] [59] shows that for any sequential ordering \(R\) of the nodes of the graph, the removal of those arcs \((i,j)\) for which \(R(i)\geq R(j)\) must eliminate all directed loops, and if all these arcs are selected as iteration variables, then, since no directed loops remain, the RISM can be evaluated in the order determined by \(R\). Hence, the chosen ordering determines the set of iteration variables. To find a minimum number of such variables, it is therefore desired to find a sequential ordering \(\mathbb{R}\) which minimizes the number of arcs \((i,j)\) for which \(R(i)\geq R(j)\). This is precisely the problem of determining a minimum feedback arc set. Details can be found in [PERSON], 1963, but also other and more efficient algorithms can be found in the literature (e.g., Ref. [60]), giving sometimes optimal but often suboptimal ordering. The general problem is an NP-hard problem [61]. Many suboptimal ordering algorithms have the great advantage of fast execution even only directly proportional to \(N\). The proposed method is analysed in detail for LTI systems in the next paragraph. The New Approach to Individual Subsystems Models: Linear Time-Invariant Systems and State-Space Representation  
The partitioning method outlined in the previous paragraph is now applied to LTI systems; Equations (3a), (3b) and (4) can be written explicitly. To distinguish vectors from matrices, we use uppercase letters for latter and lowercase letters for the former (e.g., \(x\) instead of the \(X\) symbol used in the previous paragraph):  
\[\dot{x}=Ax+Bv+Pu, \tag{9a}\] \[y=Cx,\] (9b) \[v=My, \tag{9c}\]  
where matrices \(A\), \(B\), and \(C\) are block diagonal, each block representing each subsystem in the ISM: they can also depend on time and in this case the system is not autonomous, making the discussion only more complicated. Matrix \(P\) represents the dependence of the status vector \(x\) on the external input vector \(u\). Moreover, the explicit dependence of internal outputs \(y\) on internal inputs \(v\) and external inputs \(u\) has been neglected to simplify the following calculations without loss of generality, but it can be easily reintroduced (for example for subsystems without ODEs but only algebraic input--output blocks: \(y_{i}~{}=~{}D_{ii}u_{i}+E_{ii}v_{i}\) for the i-th subsystem). Regarding external outputs, they can be chosen without loss of generality among the components of the internal output vector \(y\). Substituting the last two Equations (9b) and (9c), into the first one (9a) we obtain the following system of equations that can be solved by calculating \(x(t)\) as all the matrices and external inputs \(u\) are known  
\[\dot{x}=Kx+Pu, \tag{10}\]where the matrix \(K=A+BMC\) corresponds to the state transition matrix (precisely for LTI systems, the state transition matrix is \(e^{K\!I}\)). The system in Equation (10) is said asymptotically stable (AS) if for each \(\varepsilon>0\) there is \(\delta(\varepsilon)\) so that for \(u=0\) and a suitable norm \(\|\cdot\|\) we have  
\[\|x(0)\|<\delta\Rightarrow\left\{\|x(t)\|<\varepsilon\;\forall t>0\right. \tag{11}\]  
The system is AS if and only if the real parts of the eigenvalues \(\lambda(K)\) of matrix \(K\) are all negative: \(\lambda(K)<0\). If the sequence of discrete-time points \(t_{n}\) is equally spaced (\(t_{n}=n\Delta t\)) the discrete version of system 10 can be obtained approximating \(\dot{x}\) as follows  
\[\dot{x}\sim\frac{x_{n+1}-x_{n}}{\Delta t}, \tag{12}\]  
where \(x_{n}=x(t_{n})\). With no loss of generality, other approximations may be used to improve convergence, especially when dealing with stiff problems.